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Abstract 


A  random-choice  method  (RCM)  is  developed  for  obtaining  fairly 
practical  and  efficient  numerical  solutions  for  two-dimensional  planar 
and  axisymmetric  steady  supersonic  flows#  such  as  those  for  sharp- 
edged  planar  airfoils,  supersonic  inlets  of  aircraft  engines,  pointed 
bodies  of  revolution,  supersonic  nozzles,  and  free  jets.  This  method 
is  based  on  the  solution  of  a  Riemann  problem,  which  is  the  elemental 
solution  of  the  hyperbolic  equations  of  two-dimensional  steady  super¬ 
sonic  flows.  The  Riemann  problem  consists  of  two  waves  separated  by 
a  slip  stream,  and  each  wave  can  be  either  an  oblique  shock  wave  or  a 
Prandtl-Meyer  expansion  wave.  Advanced  techniques  are  given  for  solv¬ 
ing  the  Riemann  problem  iteratively,  handling  the  boundary  conditions 
along  body  and  free-jet  surfaces,  and  computing  only  certain  parts  of 
flow  fields  of  interest.  Many  interesting  and  practical  numerical 
solutions  are  presented  for  different  types  of  planar  and  axisymmetric 
flows,  to  demonstrate  the  applicability,  capability,  and  limitations 
of  the  RCM.  Numerical  results  are  shown  to  be  in  excellent  agreement 
with  both  known  analytical  solutions  and  results  from  the  method  of 
characteristics . 
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Appendix  A:  COMPUTER- PROGRAM  LISTING  FOR  SOLVING 
TWO-DIMENSIONAL  STEADY  SUPERSONIC 
FLOW  PROBLEMS .  A1 


v 


a 


A 

B 

C 

CP 

e 

f(W) 

g(W) 

G 

h(W,r) 

i 

j 

M 

P 

P(x,r) 

r 

rb 

R 

S 

t 

T 

u 

v 

W 

x 

Xb 


gas  sound  speed 
constant  defined  by  Eq.  46 
constant  defined  by  Eq.  46 
constant  defined  by  Eq.  46 
pressure  coefficient 

total  energy  per  unit  volume  [e  =  e  +  p(u2  +  v2)/2] 
function  of  W  in  Eqs .  1  and  2 
function  of  W  in  Eqs.  1  and  2 
vector  defined  by  Eq.  45 

function  of  W  and  radius  r  in  Eqs.  1  and  2 

integer  for  numbering  mesh  points  in  the  axial  direction 

integer  for  numbering  mesh  points  in  the  radial  direction 

flow  Mach  number 

static  pressure 

sampling  point  for  the  Riemann  problem 

radial  coordinate 

radius  of  an  axisymmetric  body 

gas  constant 

a  state  for  the  Riemann  problem 

thickness  of  a  planar  airfoil 

static  gas  temperature 

axial  component  of  the  gas  velocity 

radial  component  of  the  gas  velocity 

solution  vector  [p,  pu,  pv,  e]  from  Eqs.  1  and  2 

axial  coordinate  along  the  free-stream  flow 

locus  of  the  airfoil  surface 


Greek  symbols 


u  coefficient  for  planar  (a  =  0)  and  axisymmetric  (a  =  1)  flows 

0  angle  of  a  characteristic  line  or  a  Mach  wave  (measured 

clockwise  from  the  radial  axis) 

Y  specific-heat  ratio 

0  angle  of  flow  deflection  through  a  shock  or  rarefaction  wave 

^b  boundary— sur face  angle  measured  from  the  longitudinal  axis 

A  axial  and  radial  increments  (e.g..  Ax  and  Ar) 


Notation  (continued) 


e 

n 

K(  T| ) 

e 

«c 

p 


V(M) 

5 


n 

P 


internal  energy  per  unit  volume 
pressure  ratio  p,/pj  and  p„/pr 
function  of  q  defined  by  Eq.  16 
angle  of  a  radial  line  from  a  cone  apex 
semivertex  angle  of  a  wedge  or  cone 

Mach  angle  [p  =  sin“^(l/M)],  also  a  function  of  q  and  Ma 
from  Eq.  17 

Prandt 1-Meyer  function  in  Eq.  8 

flow  angle  n/2  -  tan'^(v/u) ,  measured  clockwise  from  the 
radial  axis  (also,  angle  of  a  slip  surface  or  slip  stream) 

ratio  of  circumference  and  diameter  of  a  circle  (3.1415927) 

static  density 

shock  angle  measured  clockwise  from  the  radial  axis 


X(q,Ma)  function  of  q  and  Ma  from  Eq.  18 

Q  value  of  the  random  variable  equidistributed  in  the 

range  -1/2  to  +1/2 


Subscripts 

a  denotes  conditions  ahead  of  an  oblique  shock  or  expansion 

wave 

b  denotes  conditions  behind  an  oblique  shock  or  expansion  wave 

h  denotes  the  value  at  the  rarefaction-wave  head 

j  index  for  numbering  the  radial  position  of  a  grid  point 

1  denotes  the  initial  state  of  the  gas  on  the  left  side  of 

the  discontinuity  in  the  Riemann  problem 

min  minimum  value 

r  denotes  the  initial  state  of  the  gas  on  the  right  side  of 

the  discontinuity  in  the  Riemann  problem 

t  denotes  the  value  at  the  rarefaction-wave  tail 

*  denotes  a  common  value  on  either  side  of  the  slip  surface 

Superscripts 

i  index  for  numbering  the  axial  position  of  a  grid  point 

overhead  bar  denotes  the  sampled  solution  for  the 
planar-flow  case 

overhead  tilde  denotes  the  first-order  sampled  solution 
in  MacCormack's  operator-splitting  scheme  (Eq.  48) 
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1.0 


INTRODUCTION 


1 . 1  Literature  Review 

Numerical  models  and  analyses  have 
simulation  of  complex  fluid-dynamic  flows, 
can  provide  an  accurate  description  of  the 
provide  a  means  of  validating  new  physical 
providing  fundamental  knowledge  about  flow 

Two-d imens iona 1 
can  be  described,  with 
vat  ion  laws.  Numerical 
following  criteria:  1) 
smooth  parts  of  the  flow, 
and  slip  streams  (or  slip 
propagated  or  transported 
should  be  computed  stably 


proven  very  powerful  for  the  study  or 
If  properly  validated,  these  models 
flow  field.  Additionally,  they  can 
and  mathematical  hypotheses,  thereby 
problems  of  interest. 


planar  and  axisymmetric  steady  supersonic  flows  of  gases 
suitable  assumptions,  by  a  hyperbolic  system  of  conser- 
solutions  of  the  hyperbolic  system  should  satisfy  the 
the  computed  solution  should  be  sufficiently  accurate  in 
2)  discontinuities  in  the  flow,  such  as  oblique  shocks 
surfaces),  should  not  only  remain  sharp  but  also  be 
at  the  correct  speed,  and  3)  such  discontinuities 
without  numerical  oscillations. 


The  method  of  characteristics  has  been  used  most  widely  up  to  the  1970s 
for  obtaining  numerical  solutions  of  steady  supersonic  flows  over  planar  and 
axisymmetr ic  bodies.  Isenberg  and  Lin* ,  Clippinger,  Giese  and  Carter^,  Ferri^, 
Shapiro'*,  and  Liepmann  and  Roshko^  found  many  important  applications  of  this 
method  for  supersonic  flows  around  thin  sharp-edged  airfoils  and  slender  pointed 
bodies  of  revolution  with  attached  shock  waves,  and  also  for  supersonic  flows 
in  aircraft-engine  inlets  and  in  many  types  of  nozzles.  In  the  application  of 
the  method  of  characteristics,  flow  properties  are  constrained  by  ordinary  dif¬ 
ferential  equations  that  apply  along  characteristic  lines  in  flow  regions  with¬ 
out  discontinuities,  so  that  solutions  can  be  obtained  on  a  characteristic  mesh. 
Special  adaptive  procedures  are  needed  to  include  and  track  both  oblique  shocks 
and  slip  streams,  and  this  is  the  primary  difficulty  in  employing  this  method. 
Each  flow  problem  normally  needs  a  new  computer  program,  because  the  computer 
logic  is  specific  to  the  problem.  As  a  result,  the  method  of  characteristics  is 
not  often  used  to  solve  flow  problems. 

The  finite-difference  method^-®  is  another  numerical  means  of  solving 
flow  problems  which  has  been  widely  used.  In  order  to  incorporate  shocks  and 
slip  streams  in  the  solutions,  special  schemes  like  artificial  and  numerical 
viscosity  have  to  be  employed.  For  example,  von  Neuman  and  Richtmyer^  were  the 
first  to  introduce  an  artificial  viscosity  term  into  the  Lagrangian  form  of  the 
gasdynamic  equations.  This  scheme,  like  those  and  others  that  followed,  auto¬ 
matically  handle  shocks  and  slip  streams,  but  they  either  introduce  undesirable 
oscillations  in  the  solution  near  these  discontinuities  or  smear  them  out  over 
several  grid  points.  Automatically  handling  discontinuities  is  so  advantageous 
that  the  finite-difference  method  is  widely  used  today  for  all  sorts  of  flow 
problems,  despite  the  fact  that  the  undesirable  smearing  and  also  incorrect 
placement  of  discontinuities  affects  the  quality  and  accuracy  of  the  solution. 

Over  the  past  two  decades  the  finite-element  method  has  been  developed 
for  the  solution  of  flow  problems*-®.  To  the  present  time,  this  method  has  been 
used  frequently  to  solve  successfully  problems  involving  incompressible  flows 
and  steady  compressible  subsonic  flows,  but  it  has  not  been  sufficiently  well 
developed  at  present  to  solve  hyperbolic  equations  with  similar  success. 

The  particle-in-cell  (PIC)  and  fluid-in-cell  (FLIC)  methods,  have  been 
applied  with  success  to  solve  a  wide  variety  of  multifluid  flow  problems  that 


1 


involve  large  slippages  and  distortions  lifce  those  that  take  place  when  a 

projectile  penetrates  a  liquid  or  metal.  The  PIC  method  is  a  combined  Eulerian 
and  Lagrangian  scheme  that  utilizes  Lagrangian  fluid  particles  to  transport 
mass,  momentum,  and  energy  through  an  Eulerian  mesh  of  cells, ^  whereas  the  FLIC 
method  uses  only  an  Eulerian  scheme  to  transport  fluid.  12-13  Although  the  use 
of  such  transported  particles  or  fluid  facilitates  in  obtaining  the  solution  of 
a  problem,  these  methods  are  not  widely  employed  because  the  predictions  contain 
aphysical  fluctuations  in  fluid  quantities  and  excessive  smearing  of  disconti¬ 
nuities,  and  they  also  normally  demand  a  large  computer  memory  and  computation 
t  ime  . 


Taylor  and  Maccoll^  and  Maccoll^  in  their  early  studies  (1933,  1937) 
solved  the  flow  around  an  infinite  light-circular  cone  moving  supersonically 
with  an  attached,  conical,  bow  shock  of  constant  strength.  For  this  special 
problem,  they  were  able  to  obtain  a  self-similar  solution  in  terms  of  ordinary 
differential  equations,  and  extensive  tables  of  flow  properties  around  cones, 
which  is  based  on  these  early  studies,  can  be  found  in  Ref.  16. 

If  a  body  in  a  supersonic  flow  is  sufficiently  slender  and  smooth,  then 
flow  deflections  or  perturbations  will  be  small.  Shocks  will  be  sufficiently 
weak  such  that  the  flow  for  practical  purposes  can  be  considered  to  be  homen- 
tropic  (i.e.,  uniform  entropy)  and  therefore  irrotat  ional  .  For  this  case  the 
most  popular  theories  introduce  a  velocity  potential  to  obtain  an  approximate 
linear  perturbation  solution.  The  perturbation  flow  solution  for  disturbance 
produced  by  the  body  is  defined  such  that  it  is  simply  superposed  on  the  uniform 
free-stream  velocity.  Such  perturbation  solutions  can  be  found  in  papers  by  von 
Karman  and  Moore-^,  Lighthilll^-19 ^  Ward^®,  Whitham^l~22 ,  an(j  Van  Dyke^,  and 
they  have  been  used  recently  by  Ritzel  and  Gottlieb^  to  obtain  the  overpressure 
signatures  from  actual  projectiles  and  by  Devan^  to  get  surface  pressures  on 
axisymmetr ic  bodies.  Although  perturbation  solutions  are  easy  to  use,  they  are 
both  approximate  and  limited  to  supersonic  flows  of  low  free-stream  Mach  number. 

A  relatively  new  method  has  been  developed  and  used  in  recent  years  for 
solving  hyperbolic  equations,  provided  that  the  elemental  Riemann  problem  has  a 
known  solution.  The  random-choice  or  Glimm's  method  has  the  highly  desirable 
capability  of  automatically  including  and  correctly  placing  discontinuities  like 
shocks  and  contact  surfaces  in  the  numerical  solutions,  without  any  artificial 
or  numerical  viscosity,  or  any  other  special  shock  capturing  techniques.  It  was 
originally  developed  for  solving  one-dimensional  unsteady  planar  flow  problems 
by  Glimm^v,  based  to  some  extent  on  earlier  work  of  Godunov^?.  Glimm's  method 
was  made  more  practical  for  solving  planar  shock-tube  and  detonation  problems 
by  Chorin^®  first  and  later  by  Co  le  1  laA^-^® ,  an(j  also  extended  by  means  of  an 
operator-splitting  technique  to  solve  cylindrical  and  spherical  flow  problems 
by  Sod^l  and  quas  i-one-d  imens  iona  1  gas  flows  by  Fok^,  Greatrix  and  Gottlieb^, 
and  Glimm,  Marshall  and  Plohr^.  Many  different  one-dimensional  unsteady-flow 
problems  have  been  solved  successfully  in  recent  years  at  UTIAS.35-48 

A  well-known  mathematical  analogy  exists  between  two-dimensional  super¬ 
sonic  steady  flows  and  one-dimensional  unsteady  flows.  The  essence  of  this 
analogy  is  that  both  problems  are  described  by  hyperbolic  systems  of  conserva¬ 
tion  laws  in  two  independent  variables,  and  each  has  a  well-known  elemental 
Riemann  problem  and  analytical  solution.  Since  the  initial  success  of  develop¬ 
ing  and  employing  the  random-choice  method  (RCM)  to  obtain  solutions  for  one¬ 
dimensional  unsteady  flow  problems,  it  was  recognized  that  similar  success  might 
also  be  achievable  for  solving  analogous  two-dimensional  steady  supersonic  flow 
problems,  provided  an  analogous  RCM  could  be  developed. 
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Marshall  and  Plohr^  were  the  first  to  develop  a  RCM  for  solving  two- 
dimensional  steady  supersonic  flow  problems  having  both  planar  and  ax i symmetric 
geometries.  At  the  same  time.  Honma,  Wada  and  Inomata®®  developed  independently 
a  similar  method.  In  both  cases,  the  Riemann  problem  was  posed  and  its  solution 
given  in  a  very  similar  manner.  However,  Marshall  and  Plohr  used  a  rectangular 
mesh  in  their  computations,  whereas  Honma,  Wada  and  Inomata  used  combined  body- 
fitted  and  rectangular  meshes.  Although  this  combined  mesh  system  might  appear 
as  beneficial  initially,  it  actually  increases  computation  time,  is  difficult  to 
use  with  a  curved  body  surface,  and  is  also  unnecessary. 

In  the  studies  by  Marshall  and  Plohr^,  and  Honma,  Wada  and  Inomata^®, 
most  of  the  effort  was  expended  on  developing  the  RCM,  and  only  a  few  simple 
exemplary  problems  were  solved  to  illustrate  feasibility.  In  the  first  study, 
numerical  solutions  were  given  for  the  supersonic  flow  over  straight  and  curved 
wedges,  and  also  a  cone  and  cylinder  combination.  In  the  other  study,  solutions 
were  given  for  the  flow  over  a  double  wedge,  a  cone,  and  a  double  cone. 


1 .2  Scope  of  the  Present  Stud^ 


The  impetus  for  the  present  study  stems  originally  from  the  experimental 
and  analytical  work  of  Ritzel  and  Gottlieb24»51-54 ,  An  exact  numerical  solution 
for  the  flow  field  around  supersonic  projectiles  was  being  sought  to  complement 
the  existing  approximate  analytical  method, ^  to  explain  and  accurately  predict 
many  features  of  the  measured  overpressure  s ignatures  .51-54  jt  appeared  that 
the  RCM  would  be  an  ideal  method,  if  it  could  be  developed  for  two-dimensional 
steady  supersonic  flows,  based  on  the  previous  success  with  the  RCM  for  solving 
one-dimensional  unsteady  flow  problems  at  UTIAS .33 » 35-48  Qur  development  of  a 
RCM  for  solving  two-dimensional  steady  supersonic  flow  problems  began  just  as  we 
became  aware  of  the  developments  by  Marshall  and  Plohr^  and  also  by  Honma,  Wada 
and  Inomata^®,  and  this  former  work  became  very  helpful  to  us. 


In  this  report  the  RCM  is  presented  and  explained  in  detail  for  both 
reference  and  completeness,  even  though  the  RCM  analysis  was  given  briefly  but 
elegantly  in  the  previous  work.  A  physical  description  is  also  presented  for 
the  elemental  Riemann  problem,  which  consists  of  two  waves  separated  by  a  slip 
stream,  and  each  wave  can  be  either  an  oblique  shock  wave  or  a  Prandt 1-Meye r 
expansion  wave.  This  analysis  is  embodied  in  a  computer  program  for  obtaining 
both  practical  and  efficient  numerical  solutions  of  two-dimensional  planar  and 
ax isymmetr ic  steady  supersonic  flows,  such  as  those  for  sharp-edged  planar  air¬ 
foils,  supersonic  inlets  of  aircraft  engines,  pointed  bodies  of  revolution, 
supersonic  nozzles,  and  free  jets.  Numerical  results  from  this  computer  program 
are  shown  to  be  in  excellent  agreement  with  both  known  analytical  solutions  and 
predictions  by  the  method  of  characteristics.  Additionally,  many  interesting 
and  practical  numerical  solutions  are  presented  for  different  types  of  planar 
and  axisymmetric  flows,  in  order  to  demonstrate  the  applicability,  capability, 
and  also  some  severe  restrictions  of  the  RCM. 


Advanced  techniques  are  presented  in  this  report  to  make  the  computer 
program  both  practical  and  efficient  for  solving  two-dimensional  planar  and 
axisymmetric  steady  supersonic  flow  problems.  Firstly,  an  efficient  parabolic- 
curve  iteration  method  is  used  in  the  solution  of  the  Riemann  problem,  making 
the  number  of  iterations  normally  three  or  less.  Secondly,  boundary  conditions 
are  handled  naturally  and  easily  without  special  coordinate  transformations.  At 
a  solid  body  the  flow  is  made  tangent  to  the  surface,  and  at  the  edge  of  a  free 
jet  the  flow  pressure  is  made  equal  to  the  ambient  pressure  outside  of  the  flow. 
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Thirdly,  an  efficient  operator-splitting  technique  is  used.  Finally,  a  scheme 
is  devised  for  numerical  computations  of  certain  problems,  so  that  only  relevant 
parts  of  the  flow  field  are  computed,  in  order  to  significantly  reduce  computa¬ 
tional  time.  For  example,  the  part  of  the  flow  field  containing  only  the  wave 
emanating  from  a  supersonic  projectile  is  computed,  and  the  free-stream  flow 
before  the  wave  and  the  slightly  varying  flow  far  behind  the  wave  is  omitted. 

The  random-choice  method  cannot  solve  all  planar  and  axisymmetric  steady 
supersonic  flow  problems  over  arbitrary  shaped  bodies.  Limitations  of  the  RCM 
are  therefore  documented  and  discussed  in  this  report.  Previous  reports  do  not 
even  mention  most  of  these  limitations. 


2.0  ANALYSIS  FOR  THE  RANDOM-CHOICE  METHOD 

2.1  General  Equations  for  Planar  and  Ax i s vmma t r ic  Flows 

The  partial  differential  equations  for  two-dimensional  steady  supersonic 
inviscid  flows  of  a  gas  with  planar  and  axisymmetric  geometry-*  can  be  written  as 

^IftW)]  +  |^Ig  (W)  ]  =  h(W,r) ,  (1) 

where 


p 

p  u 

w  = 

p  tt 

II 

p  u2  +  p 

P  V 

p  n  v 

e 

u(e  +  p) 

P  v 

P  v 

g(W)  = 

p  u  V 

p  v2  +  p 

• 

h(W,r)  =  -  f 

p  u  V 

p  V2 

v(e  +  p) 

v(e  +  p) 

The  symbols  p  and  p  are  the  static  pressure  and  density  of  the  gas,  u  and  v  are 
the  axial  (x)  and  radial  (r)  components  of  the  flow  velocity,  pu  and  pv  are  the 
axial  and  radial  components  of  momentum,  and  a  is  1  for  an  axisymmetric  flow. 

In  the  case  of  a  planar  flow,  the  axial  distance  x  becomes  the  longitudinal 
distance  x,  the  radial  distance  r  becomes  the  lateral  distance  y,  and  o  =  0  for 
this  case.  Hence,  in  the  case  of  planar  flows  the  nonhomogeneous  term  h(W,r) 
does  not  appear  and  the  equations  are  homogeneous. 

The  total  energy  per  unit  volume  e  in  Eqs .  1  and  2  is  given  by 

e  =  e  +  p(u2  +  v2)/2  (3) 

and  contains  the  internal  energy  per  unit  volume 

a  =  p/(y  -  1) ,  (4) 
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where  y  denotes  the  specific  heat  ratio  (constant  for  a  perfect  gas).  Note  that 
implicit  in  these  equations  is  the  equation  of  state  p  =  pRT  and  the  sound-speed 
relation  a^  =  yRT. 

The  four  expressions  that  are  embodied  in  Eqs .  1  and  2  are  the  continu¬ 
ity  equation,  Newton's  second  law  of  motion  of  momentum  in  the  axial  and  radial 
directions  for  an  axisymmetric  flow  (or  longitudinal  and  lateral  directions  for 
a  planar  flow),  and  the  conservation  of  energy.  This  system  of  equations,  that 
is  Eqs.  1  to  4,  can  be  solved  for  a  particular  problem  only  if  the  appropriate 
initial  and  boundary  conditions  are  specified. 


2.2  General  Riemann  Problem 


For  the  solution  of  the  equations  presented  in  the  preceding  section  by 
means  of  the  random-choice  method  (RCM),  a  Riemann  problem  must  be  established 
first  and  subsequently  solved.  The  general  Riemann  problem  is  considered  in 
this  section,  and  the  particular  solution  for  planar  and  axisymmetric  supersonic 
flows  is  given  later.  The  presentation  is  similar  to  that  of  Chorin^®. 

Consider  the  flow  field  shown  in  Fig.  1,  in  which  the  free-stream  flow 
is  aligned  approximately  along  the  axial  coordinate  x  and  the  radial  coordinate 
r  is  approximately  normal  to  the  stream  line.  Let  the  flow  field  be  subdivided 
by  means  of  a  grid  of  constant  Ax  spacing  in  the  x  direction  and  constant  Ar  in 
the  radial  direction.  Also  let  the  indicies  i  and  j  be  the  labels  for  the  grid 
corners.  Note  that  each  cell  of  width  Ax  and  height  Ar  is  further  subdivided 
into  quarters  by  the  dashed  lines,  resulting  in  alternating  grid  points  at  the 
cell  corners  and  centers  for  the  computations. 

Let  the  actual  flow  field,  which  may  be  either  smoothly  varying  in  terms 
of  its  variables  or  smoothly  varying  with  discontinuous  lines  or  surfaces,  be 
discretized  such  that  the  values  of  the  variables  of  importance  are  either  known 
initial  conditions  or  computed  ones  only  at  the  grid  corners  and  grid  centers. 

In  first-order  Riemann  initial-value  problems,  the  states  at  two  adjacent  grid 
nodes  are  not  joined  by  a  smoothly  varying  line  or  surface,  but  instead  the 
'constant*  states  are  always  joined  or  separated  by  means  of  a  discontinuity. 

The  exact  location  of  this  created  discontinuity  is  unknown,  or  the  location  of 
an  actual  one  is  unknown  if  it  existed  —  being  lost  in  the  discretization 
process.  However,  a  discontinuity  is  inserted  when  necessary  at  an  appropriate 
location  by  means  of  an  appropriate  random-sampling  process. 

When  inserting  a  discontinuity  in  the  interval  Ar,  the  only  information 
available  is  that  it  occurs  somewhere  in  this  interval.  Therefore,  one  can  only 
assume  that  it  has  an  equal  probability  of  occurring  at  each  point  or  in  each 
fraction  of  the  interval.  Hence,  a  random  variable  Q  that  is  equidistributed  in 
an  interval  from  -1/2  to  +1/2  may  be  chosen  to  get  an  appropriate  location  for 
the  discontinuity.  If  Q  =  -1/2,  for  example,  the  discontinuity  can  be  inserted 
at  r  =  jAr,  if  Q  =  0  it  can  be  inserted  at  (j+l/2)Ar,  and  if  8  =  1/2  it  can  then 
be  inserted  at  (j+l)Ar.  Note  that  this  random-sampling  process  for  inserting  or 
locating  a  discontinuity  in  the  Riemann  problem  gives  the  random-choice  method 
its  name. 

In  establishing  the  general  Riemann  problem,  let  Wj  denote  W[iAx,jAr], 
which  are  the  initial  values  or  a  past  solution  at  coordinate  point  [iAx,jAr]. 
Also,  let  Wj+i  denote  W[iAx , ( j+l)Ar] ,  which  are  also  initial  values  or  a  past 
solution  at  the  adjacent  coordinate  point  [iAx, (j+l)Ar] .  Now,  if  these  two  sets 
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of  initial  conditions  or  two  past  solutions  are  taken  as  the  initial  conditions, 
then  the  solution  Wj+J/2  farther  downstream  at  the  center  of  the  cell  or  at  the 

coordinate  point  [ ( i+l/2)Ax , ( j+l/2)Ar]  needs  to  be  determined,  to  advance  the 
solution  in  the  axial  direction.  The  initial-value  problem  which  will  give  this 
solution  can  now  be  expressed  as 


W[ iAx , r] 


+1 


if  r  >  (j  +  1/2  +  Q)Ar, 


if  r  <  (j  +  1/2  +  Q)Ar. 


(5) 


An  initial-value  problem  of  this  type  that  yields  the  solution  is  called 

a  Riemann  problem. 


The  solution  to  a  Riemann  problem  is  normally  self-similar  and  can  then 
be  expressed  in  terms  of  tlx.  When  such  a  solution  is  obtained,  it  will  extend 
outward  from  the  discontinuity  on  the  initial  line  x  =  iAx,  in  the  radial  direc¬ 
tion  r  and  downstream  in  the  direction  x,  eventually  overlapping  the  coordinate 
point  [ ( i+1/2) Ax , ( j+1/2) Ar] .  The  set  of  values  from  this  solution  that  falls 
directly  on  the  point  [ ( i+l/2)Ax , ( j +l/2)Ar]  can  then  be  assigned  as  the  solution 
mi+l/2 

Wj+i/2  to  this  grid  node. 

Once  the  solution  of  one  Riemann  problem  has  been  obtained,  it  is  not 
very  difficult  to  see  how  this  process  can  be  repeated  for  all  cells  having  the 
common  axial  distance  iAx,  in  order  to  get  the  solutions  at  all  cell  centers  at 
axial  distance  (i+l/2)Ax.  Then,  the  solution  is  advanced  column  by  column  in 
the  axial  direction,  in  order  to  obtain  the  solution  for  the  entire  flow  field. 
In  this  process,  however,  suitable  boundary  conditions  would  have  to  be  applied 
at  the  outer  edges  of  the  grid. 


It  is  worth  mentioning  here  that  the  description  of  the  random-sampling 
process  and  assignment  of  a  solution  of  the  Riemann  problem  to  the  next  grid 
node  differ  from  those  in  Chorin's  paper^S.  He  always  puts  the  discontinuity 
right  at  the  center  of  the  cell,  random  samples  the  Riemann  solution  downstream, 
and  then  assigns  this  set  of  values  to  the  next  grid  node.  Although  the  two 
methods  are  mathematically  equivalent,  the  present  method  and  description  really 
provides  more  physical  insight  and  understanding. 


2 . 3  Solution  of  the  Riemann  Problem  for  Steady  Supersonic  Planar  Flows 

The  solution  of  the  Riemann  problem  cannot  always  be  obtained  in  simple 
analytical  form  for  complex  sets  of  partial  differential  equations.  This  is, 
unfortunately,  the  case  for  steady  supersonic  axisymmetric  flows,  because  the 
inhomogeneous  term  h(W,r)  unduly  complicates  the  equations.  When  this  term  is 
absent,  as  in  the  specific  case  of  steady  supersonic  planar  flows,  an  analytical 
solution  can  then  be  obtained.  These  solutions  for  planar  flows  are  given  first 
in  this  section,  and  the  operator-splitting  technique  which  can  be  applied  to 
modify  the  planar  solutions  to  yield  those  for  axisymmetric  flows  is  introduced 
later . 


The  Riemann  problem  for  the  hyperbolic  system  of  equations  given  by  Eqs . 
1  and  4  is  an  initial  value  problem  for  which  the  initial  data  is  known  at  an 
axial  distance  iAx  (Fig.  1),  as  mentioned  earlier.  The  initial  data  consist  of 
two  states,  on  the  left  (at  the  larger  radius),  and  Sr  on  the  right  (at  the 
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smaller  radius).  These  two  states  are  separated  by  an  initial  discontinuity  at 
some  radial  distance  r  (see  Eq.  5),  and  it  breaks  into  left  and  right  running 
waves  that  are  swept  downstream  by  the  oncoming  supersonic  flow  and  separated  by 
a  slip  surface  or  slip  stream,  as  shown  in  Fig.  2. 

In  the  Riemann  problem  for  steady  supersonic  flows,  the  left  and  right 
running  waves  can  be  either  an  elemental  oblique  shock  wave  or  Prandtl-Meyer 
rarefaction  wave. 3  These  two  elemental  waves  result  in  four  different  wave 
patterns  and  solutions  (see  Fig.  2).  The  left  and  right  running  waves  are 
separated  by  a  slip  stream,  across  which  the  pressure  and  flow  angle  are  con¬ 
tinuous  but  the  density,  temperature,  flow  velocity,  and  flow  Mach  number  are 
normally  discontinuous.  For  the  case  of  steady  supersonic  planar  flows  of  a 
polytropic  gas,  the  partial  differential  equations  for  the  Riemann  problem  can 
be  reduced  to  analytical  equations,  and  these  will  soon  be  presented.  For  the 
case  of  a  polytropic  gas,  only  four  variables  are  needed  to  fully  describe  a 
particular  state,  whether  this  state  is  ahead  of  a  shock  or  rarefaction  wave 
('a'  in  Fig.  2),  behind  a  shock  or  rarefaction  wave  ('b'),  or  inside  a  rare¬ 
faction  wave.  As  in  past  studies49-50f  the  four  most  convenient  variables  are 
the  static  pressure  p,  static  density  p,  flow  Mach  number  M  =  ( {u2  +  v^J/a2)1^, 
and  flow  angle  £  =  n/2  -  tan-1(v/u) . 

Consider  the  case  of  a  rarefaction  wave  first.  In  this  case  the  partial 
differential  equations  reduce  to  ordinary  differential  equations  that  then  apply 
along  characteristic  lines,  which  can  be  further  integrated  for  planar  flows  to 
obtain  the  Prandtl-Meyer  solution. ^  The  centered  fan  of  characteristic  lines 
of  a  rarefaction  wave  is  depicted  in  Fig.  3,  for  rarefaction  waves  on  the  left 
and  right  sides  of  the  stream  line.  In  this  diagram  three  different  angles  are 
depicted.  The  Mach  angle  p  or  the  angle  of  a  characteristic  line  with  respect 
to  the  stream  line  is  given  by  sin“1(l/M) .  The  angle  of  the  flow  £  is  measured 
from  the  radial  axis  to  the  stream  line.  Therefore,  the  angle  of  a  particular 
characteristic  line,  measured  from  the  radial  axis,  is  given  by 

P  =  $  +  sin-1(l/M) ,  (6) 

where  the  minus  sign  is  used  for  the  left  rarefaction  wave  and  the  positive  sign 
is  used  for  the  right  rarefaction  wave.  The  state  of  the  gas  is  constant  along 
one  particular  characteristic  line  having  an  angle  p.  The  flow  angle  (  and  Mach 
number  M  for  this  characteristic  line  within  the  fan  of  characteristics  can  be 
related  to  the  initial  flow  properties  (with  subscript  a)  by  the  equation 

i  5  =  *  5a  +  v(M>  -  v(Ma),  (7) 

where 

V(M)  =  [T-^l]1/2tan_1[[l-^-^]1/2[M2  -  1] 1/2]  -  tan_1[M2  -  I]1/2  (8) 

is  the  well-known  Prandtl-Meyer  function^.  The  static  pressure  and  density  for 
this  characteristic  line  be  related  to  the  flow  Mach  number  by  means  of 

p/pa  =  [(1+  3LyAmJ)/(1+  ]  (9) 

and 

P  =  Pa<p/Pa>1/Y-  (10) 
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Equation  9  can  also  be  rewritten  in  the  following  form 


M 


[ 


p  Pa 
Pa  r 


(1  + 


(11) 


which  will  be  important  for  later  work.  Because  the  pressure  decreases  from  one 
characteristic  line  to  the  next  through  a  rarefaction  wave,  from  its  head  to  its 
tail,  p  <_  pa.  Also,  p  <_  pa,  and  M  >.  Ma . 


Now  consider  the  opposite  case  of  an  oblique  shock  wave.  In  this  case 
the  partial  differential  equations  also  yield  analytical  results  for  the  sudden 
changes  in  flow  properties  across  the  shock. ®  From  the  diagrams  given  in  Fig.  4 
it  should  be  obvious  that 


1  5b  =  1  5a  +  [n/2  -  os  -  6]  -  tn/2  -  o8],  (12) 

where  (a  and  4b  ace  tlie  stream  line  angles,  measured  from  the  radial  axis,  ahead 
of  and  behind  the  oblique  shock  wave,  os  is  the  oblique  shock  angle  measured 
from  the  radial  axis,  and  6  is  the  flow  deflection  angle.  Also,  the  positive 
and  negative  signs  go  with  the  shock  wave  on  the  left  and  right  of  the  stream 
line,  respectively.  Based  on  expressions  for  oblique  shock  waves,  the  previous 
equation  can  be  expressed  in  a  convenient  form  as 

-  <b  -  *  'f'V'b  -u]1"] 

-  [*  *  -  nf2] 

for  later  use.  Furthermore,  the  flow  Mach  number  behind  the  oblique  shock 
can  be  related  to  the  flow  Mach  number  Ma  in  front  of  the  shock  by 


«b  ■  -  1  if2-  <i4> 

where  the  density  p j,  behind  the  shock  is  given  in  terms  of  the  pressure  ratio 
Pb/Pa  as 

p.  =  p  l2y  +  (y+lMp./p  “  1)3/  [2y  +  (y-l)(p  /p  ~  1)1*  (15) 

d  a  D  a  o  a 

where  this  last  equation  is  a  more  recognizable  Rankine-Hugoniot  relation.  Note 
that  the  pressure  increases  from  the  front  to  the  rear  of  an  oblique  shock  wave 
and  therefore  pj,  >.  pa.  Furthermore,  pj,  >_  pa  and  Mj,  <.  Ma. 

The  equations  above  for  oblique  shock  waves  contain  two  solutions  —  a 
weak  shock  and  a  strong  shock  solution.  Only  the  weak-shock  solution  with  the 
smallest  entropy  increase  should  be  used,  because  it  is  the  only  one  observed 
in  actual  flow  fields. 

In  the  case  of  slip  streams  or  surfaces  that  separate  different  shock 
and  rarefaction  waves,  as  shown  in  Figs.  2  to  4,  the  pressure  p^  and  flow  angle 
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5b  are  both  continuous  across  such  streams  or  surfaces.  This  is  indicated  in 
Fig.  5.  However,  the  density  p^,  and  flow  Mach  number  are  both  discontinuous. 

The  two  previous  analytical  solutions  for  the  elemental  shock  and  rare¬ 
faction  waves,  and  the  conditions  that  both  the  flow  angle  and  pressure  must  be 
continuous  across  the  slip  line,  can  now  be  combined  or  patched  together  to  give 
four  individual  solutions  for  the  four  different  wave  patterns  shown  in  Fig.  2. 
The  type  of  wave  pattern  and  the  strengths  of  the  individual  waves  depend  on  the 
initial  conditions,  that  is  the  initial  left  state  Sj  (pj»  pj,  ,  and  £l)  and 
the  initial  right  state  Sr  (pr,  pr,  Mt,  and  5r) .  Regardless  of  the  initial  con¬ 
ditions,  the  four  solutions  can  be  combined  into  one  algorithm.  In  order  to  do 
this  in  an  elegant  way,  the  following  functions  involving  the  previous  expres¬ 
sions  are  introduced 


I  n  if  n  <  1 

K(n)  =  <  (16) 

I  [2y  +  (y+1) (q-1) 1/ [2y  +  (y-l)(q-l)]  if  n  >  1 


p(n.Ma) 


r  2*(n) 

L  (y-l)n 


t1 


(17) 


v[p(n.Ma)]  -  v(Ma) 


if  n  <  l 


X(t1  ,Ma)  = 


iin-ifti  +  -^[l/n  -  1]]  p-1(n.Ma) ] 

-  sin-1[[l  +  -  11]  Ma1  ]  if  1  >  1 


In  these  expressions  the  symbol  r]  is  the  pressure  ratio  across  either  an  oblique 
shock  or  expansion  wave  (i.e.,  p^/pa)  .  These  three  functions  can  be  applied 
across  the  left  and  right  waves  of  the  Riemann  problem  (see  Fig.  5),  regardless 
of  whether  each  is  a  shock  or  rarefaction  wave.  If  this  is  done  to  connect  the 
flow  angles  ahead  and  behind  the  left  and  right  waves. 


?*  =  ?j  +  X(p./p1.M1)  (19) 

and 

"  ?*  =  -  tr  +  X(p*/pr.Mr)  (20) 

are  obtained  in  terms  of  the  common  flow  angle  and  pressure  p#  on  each  side 
of  the  slip  stream.  These  two  equations  can  now  be  added  to  eliminate  £*  and 
thereby  obtain  the  implicit  equation 


?r  "  =  X(P./P1.II1)  +  X(p*/pr,Mr)  (21) 

for  the  common  pressure  pa.  This  is  the  oasic  equation  required  to  determine 
the  common  pressure  for  each  of  the  four  different  wave  patterns.  Then  the  flow 
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angle  5*  can  be  obtained  from  either  Eq.  19  or  20,  and  the  flow  Mach  numbers  M 
on  either  side  of  the  slip  stream  can  then  be  obtained  from  of  Eq.  17,  because 
p(q,Ma)  is  actually  the  Mach  number  behind  either  shock  or  rarefaction  wave. 
Finally,  the  densities  p  follow  directly  from  Eq.  16,  because  <(q)  is  actually 
the  density  ratio  for  a  shock  or  rarefaction  wave.  This  completes  the  basic 
solution  procedure  of  the  Riemann  problem. 

The  iterative  procedure  of  solving  Eq.  21  for  p#  that  is  used  in  this 
report  is  the  parabolic-curve  method^®,  which  is  a  common  numerical  technique  in 
China.  In  this  method,  a  function  G  =  £r  -  £l  -  X(pa/pi.Mj)  -  X(pa/pr,Mr)  is 
defined  and  initially  calculated  three  times  for  three  different  initial  values 
or  guesses  of  p0,  giving  (Gj.p^i),  (G2'P*2^»  an<*  A  parabolic  curve  is 

then  fitted  to  these  three  sets  of  values  of  G  and  p^,  and  the  value  of  p^  when 
G  equals  zero  is  obtained.  Then,  this  new  value  p^4  and  its  associated  value  G4 , 
and  the  two  previous  sets  of  values  of  G  and  p0,  that  is  (G2,pa2^  and  (®3'P*3), 
can  be  used  to  produce  a  new  parabolic  curve  and  estimate  of  p0.  This  procedure 
can  then  be  repeated  until  the  solution  of  p^  is  obtained  to  the  desired  degree 
of  accuracy  (e.g.,  error  less  than  10“^  percent). 

The  first  guess  of  p#  for  the  computer  code  which  was  used  in  this  study 
was  given  by  (pj  +  pr)/2.  For  the  parabolic-curve  method,  which  requires  three 
initial  guesses,  the  second  one  was  frequently  0.98p^  and  the  third  was  1.02p#. 

In  the  parabolic-curve  method,  there  are  always  two  roots  for  p^.  By 
using  the  smallest  root  for  p^,  the  weak  shock  solution  is  chosen  automatically 
over  the  strong  shock  solution.  This  is  one  notable  advantage  of  this  method. 
Another  advantage  is  that  convergence  to  the  root  is  accelerated  by  the  use  of 
a  parabola  instead  of  a  straight  line,  and  it  would  be  virtually  just  as  fast 
as  Newton's  method.  (In  Newton's  method  additional  equations  would  be  needed 
to  find  the  derivative  of  G.)  Also,  the  secant  method  would  likely  be  slightly 
slower  in  convergence  than  either  the  psrabolic-curve  or  Newton's  method. 


2 .4  Random-Sampling  Procedure 

The  solution  of  the  initial-value  or  Riemann  problem,  given  the  left  and 
right  states  Sj  and  Sr,  has  been  determined  in  the  previous  section.  This  solu¬ 
tion  includes  both  the  type  of  wave  pattern  and  the  strengths  of  the  particular 
waves.  The  task  now  is  to  use  this  solution  to  obtain  the  particular  solution 
at  the  next  grid  node,  spaced  downstream  by  a  distance  (l/2)Ax  and  midway  be¬ 
tween  the  original  two  grid  nodes  with  states  Sj  and  Sr.  This  is  accomplished 
with  a  random-sampling  procedure  that  will  now  be  described. 

In  describing  the  general  Riemann  problem  in  section  2.2,  it  was  briefly 
mentioned  that  the  initial  discontinuity  separating  constant  states  Sj  and  Sr 
was  specified  by  means  of  selecting  a  random  number  from  a  uniform  distribution 
over  -1/2  to  +1/2.  Hence,  the  particular  wave  pattern  from  the  Riemann  solution 
would  have  its  original  discontinuity  located  at  some  randomly  specified  point 
between  the  original  nodes  (circles  in  Fig.  5)  and  then  branch  outward  in  the 
downstream  direction,  as  shown  in  Fig.  5.  Then,  depending  on  which  part  of  the 
flow  pattern  overlaps  the  next  grid  node  downstream  (the  cross  in  Fig.  5),  the 
solution  at  this  point  in  the  flow  field  would  be  assigned  to  this  grid  node. 

For  example,  if  this  downstream  grid  node  was  between  the  slip  stream  and  the 
right  wave  (see  the  figure),  then  the  solution  of  this  recently  computed  state 
would  be  assigned  to  this  node.  If  this  downstream  node  was  to  the  right  or 
ahead  of  the  right  wave,  then  state  Sr  would  be  assigned  to  this  grid  node. 
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The  procedure  outlined  above  could  be  programmed  easily  for  a  computer. 
However,  it  is  easier  and  customary  to  program  an  alternate  but  fully  equivalent 
procedure.  In  this  procedure  the  initial  discontinuity  of  the  wave  pattern  is 
always  fixed  at  the  midpoint  between  the  initial  two  grid  nodes,  which  are  shown 
as  circles  and  having  states  Sj  and  Sr  in  Fig.  6.  Then  the  wave  pattern  is  sam¬ 
pled  randomly  downstream,  along  a  line  at  a  fixed  distance  x+(l/2)Ax  that  passes 
through  the  downstream  grid  node  (the  cross).  This  sampling  point  is  labelled  P 
in  the  figure.  If  the  random  number  Q  was  equal  to  -1/2,  then  point  P  would  be 
located  at  the  corner  labelled  a.  If  Q  was  equal  to  0,  then  the  point  P  would 
be  coincident  with  the  downstream  grid  point,  and  if  0  was  equal  to  +1/2,  then 
point  P  would  be  at  the  top  corner  labelled  b.  The  solution  for  point  P  from 
the  Riemann  problem,  denoted  as  state  S  or  in  terms  of  variables  as  p,  p,  M,  and 
£,  is  then  assigned  to  the  downstream  grid  point. 

In  the  random-sampling  procedure  there  are  four  basic  cases  that  must  be 
considered.  These  cases  are  summarized  below  and  then  dealt  with  in  detail. 

1.  The  sampling  point  P[x+(l/2)Ax,r+QAr]  lies  to  the  left  of  the 
slip  stream,  which  has  a  slope  dr/dx  =  cot(£#)  in  the  (x.r)-plane, 
and  the  left  wave  is  a  shock  wave,  that  is, 

QAr  2.  (1/2)  (Ax)cot(^)  and  p,  >  p., 

2.  The  sampling  point  P[x+(l/2)Ax,r+QAr]  lies  to  the  left  of  the 
slip  stream,  which  has  a  slope  cot(!;,)  in  the  (x.r)-plane,  and 
the  left  wave  is  a  rarefaction  wave,  that  is, 

QAr  2.  ( 1/2)  (Ax)cot  (£*)  and  p#  <.  p^. 

3.  The  sampling  point  P[x+(l/2)Ax,r+QAr]  lies  to  the  right  of  the 
slip  stream,  which  has  a  slope  cot(£#)  in  the  (x,r)-plane,  and 
the  right  wave  is  a  shock  wave,  that  is, 

QAr  <  ( 1/2) (Ax)cot (£^)  and  p#  >  p^. 

4.  The  sampling  point  P[x+(l/2)Ax,r+QAr)  lies  to  the  right  of  the 
slip  stream,  which  has  a  slope  cot(£#)  in  the  (x,r)-plane,  and 
the  right  wave  is  a  rarefaction  wave,  that  is, 

QAr  <  (1/2) (Ax)cot (£#)  and  p^  <.Pf. 

Case  1:  QAr  >.  (1/2)  (Ax) cot  ( )  and  p^  >  p^ 

In  this  case  the  shock-wave  angle  os ,  measured  from  the  radial  axis  to 
the  shock  wave  (see  Fig.  4),  is  given  by 

°s  =  $1  "  8in_1[^~  t1  +  1_2^Cp*/pl  "  1]^  l-  (22> 

If  the  sampling  point  P  lies  to  the  left  or  in  front  of  the  shock  wave,  that  is, 
if  QAr  >  (1/2) (Ax)cot(os) ,  then  the  values  for  thesampling  point  P  are  those  of 
the  left  state  Sj ,  or  p  =  pj,  p  =  pj,  M  =  Mj,  and  £  =  5l«  I®  this  subcase  these 
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values  are  then  assigned  to  the  downstream  grid  point  (shown  as  a  cross).  If 
the  sampling  point  P  lies  to  the  right  of  or  behind  the  shock  wave,  that  is, 
QAr  (1/2) (Ax)cot(os) ,  then  the  sampling-point  values  for  point  P  are  p  =  p#, 
P  =  Pi<<P*/Pl)»  M  =  p(p*/pl»Ml)<  and  5  =  5*-  1“  this  subcase  these  values  are 

then  assigned  to  the  downstream  grid  point. 


Case  2:  QAr  >_  (1/2)  (Ax)cot(£#)  and  p#  <_  p^ 

The  rarefaction  wave  is  bounded  on  the  left  by  its  head  and  on  the  right 
by  its  tail,  and  the  equations  for  these  characteristics  are  given  by 

0h  =  -  sin-id/Mj),  (23) 

=  S*  sin-1[p-1(p#/p1,M1)] ,  (24) 

where  pjj  and  3t  are  the  angles  of  the  head  and  tail  characteristics,  measured 
from  the  radial  axis  (see  Fig.  3).  If  the  sampling  point  P  lies  to  the  left  of 
or  in  front  of  the  rarefaction-wave  head,  that  is,  QAr  >.  (1/2)  (Ax)cot(3]1) ,  then 
the  sampling  point  values  for  P  are  simply  those  of  the  left  state  or  p  =  pj, 
p  =  pi,  M  =  Mj,  and  5  =  Si •  If  the  sampling  point  lies  to  right  of  or  behind 
the  rarefaction-wave  tail,  that  is,  QAr  <_  (1/2) (Ax)cot(3^) ,  then  the  values  for 
sampling  point  P  are  p  =  p^,  p  =  Pi<(p*/pi).  M  =  ptp^/pi.Mj),  and  £  =  £*• 

Finally,  if  the  sampling  point  lies  inside  the  characteristic  fan  of  the 
rarefaction  wave,  then  the  values  for  the  sampling  point  are  more  difficult  to 
obtain,  because  the  solution  is  iterative.  From  the  angle  of  the  characteristic 
line  passing  through  the  sampling  point  P,  the  equation 

l  -  s in~I (1/M)  =  n/2  -  tan“I(2QAr/Ax)  (25) 


relates  the  two  unknowns  £  and  M.  The  Prandtl-Meyer  equation, 

i  =  +  V(M)  -  V(MX),  (26) 

with  the  same  two  unknowns  can  be  used  to  eliminate  yielding  the  implicit 
equation  in  M 


V(M)  -  sin_1(l/M)  =  -  5i  +  V(MX)  +  n/2  -  tan_1(2QAr/Ax) .  (27) 


This  equation  is  solved  for  M  by  employing  the  parabolic-curve  method  that  was 
described  earlier.  Newton's  method  could  also  be  used,  and  it  would  probably 
give  a  faster  convergence.  Finally,  £  is  obtained  from  either  of  Eqs .  25  or  26, 
and  the  expressions 


P 


P1 


2  +  (r  -  dm2  -y/(r-l) 
2  +  (y  -  1)M2  J 


(28) 


P  =  Pl  <(P /?!>» 


(29) 


give  the  final  values  for  p  and  p. 


Case  3:  QAr  <  (1/2) (Ax)cot(£,)  and  pm  >  pf 

Case  3  is  a  mirror  image  of  case  1.  The  same  logic  applies  but  there 
are  some  changes  that  occur  in  the  equations,  and  these  are  mainly  subscripts 
and  mathematical  symbols. 


In  this 

the  shock  wave 

case 

(see 

the 
Fig  . 

shock-wave  angle 
4)  ,  is  g  iven  by 

as,  measured 

from  the  radial 

axis  to 

*s  = 

5r 

+ 

sin_1[li;  E1  + 

Y  +  1,  , 

2y  [P*/Pr 

-11]  ]. 

(30) 

If  the  sampling  point  P  lies  to  the  right  (in  front  of  the  shock  wave),  that  is, 
f  if  QAr  <  ( 1/2) (Ax)cot (os ) ,  then  the  values  for  the  sampling  point  P  are  those  of 

the  right  state  Sr,  or  p  =  pr,  p  =  pr,  M  =  Mr,  and  5  =  Sr-  1°  this  subcase  these 
values  are  then  assigned  to  the  downstream  grid  point  (shown  as  a  cross).  If 
the  sampling  point  P  lies  to  the  left  of  or  behind  the  shock  wave,  that  is, 

QAr  >.  (1/2)  (Ax)cot(ors)  ,  then  the  sampling-point  values  for  point  P  are  p  =  p#, 
p  =  prK(p*/pr) .  M  =  p(p^/pr,Mr),  and  §  =  In  this  subcase  these  values  are 

then  assigned  to  the  downstream  grid  point. 


Case  4:  QAr  <  (1/2)  (Ax)cot(^*)  and  p^  <_  p^ 


Case  4  is  a  mirror  image  of  case  2.  The  same  logic  applies  but  some 
changes  occur  in  the  equations,  such  as  subscripts  and  mathematics  symbols. 

The  rarefaction  wave  is  bounded  on  the  right  by  its  head  and  on  the  left 
by  its  tail,  and  the  equations  for  these  characteristics  r ’•e  given  by 

Ph  =  6r  +  sin~1(l/Mr)  (31) 

and 

Pt  =  +  sin-1 [p-1 (P*/pf ,Mf) ] ,  (32) 

where  Pj,  and  are  the  angles  of  the  head  and  tail  characteristics,  measured 

from  the  radial  axis  (see  Fig.  3).  If  the  sampling  point  P  lies  to  the  right  of 

or  in  front  of  the  rarefaction-wave  head,  that  is,  QAr  <_  (1/2) (Ax)cot (0j,) ,  then 
the  sampling  point  values  for  P  are  those  of  the  right  state  Sr,  or  p  =  pr, 

P  =  Pr*  M  =  Mr,  and  $  =  !jr.  If  the  sampling  point  lies  to  left  of  or  behind 
the  rarefaction-wave  tail,  that  is,  QAr  >_  ( 1/2 )  ( Ax)  cot  ( (3t ) ,  then  the  values  for 
sampling  point  P  are  p  =  p^,  p  =  pri<(p,/pr),  M  =  p(p,/pr,Mr),  and  5  =  5,. 

Finally,  if  the  sampling  point  lies  inside  the  characteristic  fan  of  the 
rarefaction  wave,  then  the  values  for  the  sampling  point  are  more  difficult  to 
obtain,  because  the  solution  is  iterative.  From  the  angle  of  the  characteristic 
line  passing  through  the  sampling  noint  P,  the  equation 

5  +  sin-1(l/M)  =  n/2  -  tan-1(2QAr/Ax)  (33) 

relates  the  two  unknowns  ^  and  M.  The  Prandtl-Meyer  equation  exressed  in  the 
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following  form. 


5  =  5r  -  V(M)  +  V(Mr),  (34) 

with  the  same  two  unknowns,  can  be  used  to  eliminate  yielding  the  implicit 
equation  in  M, 

v(M)  -  sin-1(l/M)  =  4r  +  V(Mr)  -  n/2  +  tan-1(20Ar/Ax) .  (35) 

This  equation  is  solved  for  M_by  employing  the  parabolic-curve  method  that  was 
described  earlier.  Finally,  \  is  obtained  from  either  of  Eqs .  33  or  34,  and 
the  expressions 

r  2  +  (y  -  1)M2  (y-D 

p  =  p  IT  I  •  (36) 

L  2  +  (y  -  1)M2  J 

p  =  Pr  <(p/pf).  (37) 


give  the  final  values  for  p  and  p. 

In  the  random-choice  method,  it  is  natural  to  think  of  using  a  new  value 
of  the  random  number  Q  for  each  cell,  that  is  for  each  combination  of  i  and  j. 
However,  the  practical  effect  of  such  a  choice  with  finite  spacings  of  Ax  and  Ar 
is  disastrous,  except  for  flow-field  data  that  is  nearly  constant,  as  originally 
pointed  out  by  Chorin.2*  If  this  is  in  fact  done  there  is  a  finite  probability 
that  a  given  state  will  propagate  to  both  the  left  and  right,  thereby  creating 
a  spurious  constant  state.  The  numerical  results  will  become  less  accurate  and 
more  jagged  (as  if  the  results  contain  numerical  noise).  In  the  random-choice 
method  used  here,  a  new  random  number  Q  is  chosen  only  once  per  new  level  in 
the  axial  direction,  that  is  once  for  each  new  step  Ax.  This  random  number  is 
then  used  for  all  cells  of  this  column.  This  is  now  a  common  practice  for  the 
random-choice  method. 


2 . 5  Random-Number  Algorithm 

The  type  of  random-number  algorithm  plays  a  significant  role  in  both  the 
behavior  of  the  solution  and  the  quality  of  numerical  results  from  the  random- 
choice  method.  Better  random  numbers  produce  numerical  results  that  are  more 
accurate  (e.g.,  more  accurate  placement  of  the  position  of  shocks  in  the  flow 
field),  and  these  results  also  have  less  jaggedness  or  numerical  noise .28-30,47 
It  is  now  known  that  the  best  random  numbers  for  the  random-choice  method  are 
.  aes  that  are  actually  nonrandom  but  equidistributed ,  or  ones  that  become  equi- 
distributed  relatively  quickly  or  alternatively  tend  as  fast  as  possible  to 
approximate  equipart it  ion  in  the  range  from  -1/2  to  +1/2.28  In  this  regard  it 
appears  at  the  present  time  that  the  best  random-number  algorithm  for  use  with 
the  random-choice  method  is  that  due  to  Van  der  Corput ,29-30,47  This  random- 
number  algorithm  appears  in  Van  der  Corput's  original  work3**,  is  given  in  the 
book  by  Hammersley  and  Handscomb32,  and  is  also  discussed  by  Colella2^-3**. 

This  algorithm  is  also  presented  here  because  it  is  rather  short  and  also  for 
the  sake  of  completeness. 
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For  the  Van  der  Corput  algorithm,  suppose  that  the  natural  numbers  are 
expressed  in  the  scale  of  notation  with  a  radix  equal  to  2,  so  that 


n 


(38) 


This  is 
where  i 
bers 
reverse 


also  the  binary  expansion  of  the  number  sequence  n  =  1,  2,  3,  .  .  ,  etc, 

is  a  binary  number  that  equals  only  0  and  1.  Then,  a  set  of  random  nan- 
can  be  obtained  by  simply  writing  the  digits  of  these  numbers  in  their 
order,  preceded  by  a  decimal  point.  This  gives  the  following  series 


Q 

n 


2-<k+1) 


(39) 


The  method  of  getting  the  sequence  of  random  numbers  from  the  above  equations 
might  seem  to  be  a  little  unclear.  In  order  to  clearly  illustrate  the  manner 
in  which  this  sequence  is  constructed,  the  first  few  elements  of  the  various 
sequences  involved  are  written  down  for  convenience  in  the  table  given  below. 


n  (decimal) 

ifc  (binary) 

Qn  (binary) 

Q,!  (decimal) 

1 

1 

0.1 

0.5000 

2 

10 

0.01 

0.2500 

3 

11 

0.11 

0.7500 

4 

100 

0.001 

0.1250 

5 

101 

0.101 

0.6250 

6 

110 

0.011 

0.3750 

7 

111 

0.111 

0.8750 

8 

1000 

0.0001 

0.0625 

9 

1001 

0.1001 

0.5625 

10 

1010 

0.0101 

0.3125 

The  decimal  numbers  n  are  first  changed  into  the  binary  numbers  i^ .  Then  these 
binary  numbers  are  reversed  and  a  point  or  decimal  is  put  in  the  front  to  get  Gn 
as  binary  numbers.  Finally,  these  binary  numbers  are  simply  converted  back  into 
decimal  numbers,  yielding  Qn  in  decimal  notation,  and  covering  the  entire  range 
from  0  to  1 . 

It  can  be  noted  that  Qn  is  less  than  1/2  if  n  is  even,  and  Qn  is  greater 
than  1/2  if  n  is  odd,  but  ttjj  always  covers  the  range  from  0  to  1 .  Furthermore, 
we  have  (k/4)  <.  <  ((k+l}/4)  if  n  =  mod  (4),  with  k  =  1,  2,  and  4,  where 

JO  =  0*  jl  =  2,  J2  =  1,  and  j3  =  3.  Finally,  note  that  this  Van  der  Corput 
random-number  algorithm,  like  all  of  the  others  with  radix  different  than  2,  is 
equidistributed  . 
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2 .6  Accuracy  of  the  Flow-Field  Computations 

The  accuracy  of  computations  of  the  flow  field  depends  mainly  on  the 
size  of  the  grid.  The  constant  grid  spacing  Ar  in  the  radial  direction  is 
specified,  independent  of  the  problem.  This  spacing  depends  on  the  particular 
problem  being  solved,  the  quality  of  numerical  results  desired,  and  the  time  and 
cost  of  the  computations.  A  suitably  small  grid  spacing  in  the  radial  direction 
will  ensure  almost  any  accuracy  desired,  provided  that  computer  roundoff  errors 
are  not  encountered. 

The  grid  spacing  Ai  in  the  axial  direction  should  not  be  specified  inde¬ 
pendently  or  arbitrarily,  because  this  spacing  can  affect  the  numerical  accuracy 
of  the  results.  If  Ax  is  very  small,  numerical  accuracy  is  maintained  but  the 
number  of  computations  and  resulting  computational  cost  are  high.  Hence,  given 
the  spacing  Ar,  the  other  spacing  Ax  should  be  made  as  large  as  possible,  with¬ 
out  making  the  numerical  solution  unstable  or  inaccurate.  In  the  random-choice 
method,  which  is  an  explicit  method,  there  is  no  essential  criterion  for  numeri¬ 
cal  stability.  In  other  words,  the  solution  for  any  length  A.  is  stable.  How¬ 
ever,  if  Ax  is  too  large  the  numerical  results  will  be  inaccurate,  because  the 
effects  of  overlapping  of  adjacent  wave  patterns  is  ignored.  Experience  shows 
that  numerical  accuracy  can  be  maintained  in  the  solution  of  nonstationary  one¬ 
dimensional  flows  if  the  time  step  is  limited  by  or  specified  according  to  the 
Cou ran t -Fr iedr ic hs-Lewy  cr iterion  .28,48  This  same  criterion  is  used  here  for 
two-dimensional  steady  supersonic  flows. 

The  Courant-Friedr ichs-Lewy  criterion  is  a  limitation  on  the  size  of  Ax. 
For  this  condition  to  be  satisfied,  the  wave  pattern  cannot  spread  out  radially 
by  more  than  a  distance  of  Ar  over  the  distance  Ax.  This  means  that  adjacent 
wave  patterns  should  not  interfere  or  overlap,  otherwise  the  numerical  accuracy 
will  decrease.  This  criterion  is  expressed  as 

Ax  <.  Ar  |tan(£  +  p)lmin,  (40) 

where  £  is  the  angle  of  a  stream  line  measured  from  the  radial  axis  and  p  is  the 
Mach  angle  measured  between  the  stream  line  and  a  Mach  wave.  The  negative  sign 
is  employed  for  the  case  of  a  left  Mach  wave  and  ihe  positive  sign  is  used  with 
a  right  Mach  wave.  In  the  this  study,  the  largest  Ax  is  normally  used,  that  is. 

Ax  =  Ar  ltan(£  +  p)lmin.  (41) 

However,  whenever  Ax  is  larger  than  Ar,  which  occurs  for  highly  supersonic  flows 
because  the  wedge-shaped  pattern  is  not  widely  spread.  Ax  is  further  restricted 
by  being  set  equal  to  Ar.  This  helps  produce  more  smoothly  varying  numerical 
results  or  profiles  in  the  axial  direction. 


2 . 7  Boundary  Conditions 

Boundary  conditions  must  normally  be  specified  at  the  outer  extremities 
of  the  flow  field.  These  are  normally  needed  for  the  outermost  grid  points  in 
the  radial  direction  and  sometimes  near  or  at  the  origin.  One  or  both  of  these 
boundaries  can  be  the  direct  result  of  a  solid  body  or  a  rigid  surface  in  the 
flow,  which  are  normally  approximated  as  local  straight-line  segments  if  the 
surface  is  curved.  One  of  these  boundaries  can  also  be  an  arbitrary  edge  of 
the  flow  field,  across  which  waves  can  move  freely  without  reflection  as  they 
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leave  the  flow  field.  Finally,  the  boundary  aay  also  be  that  of  a  jet  of  gas 
or  free  jet  passing  through  another  gas  that  is  essentially  stationary  and  has 
a  basically  uniform  pressure.  All  of  these  boundary  conditions  are  considered 
here  in . 


Firstly,  consider  the  continous  boundary  condition.  In  this  simple  case 
the  flow  field  is  ended  at  a  radial  distance  that  is  specified  for  convenience. 
Any  compression,  shock,  or  expansion  wave  which  impinges  on  this  boundary  must 
be  allowed  to  pass  out  of  the  flow  field  without  any  significant  reflection  or 
influence  on  the  flow  field,  just  as  if  the  flow  field  continued  indefinitely  in 
the  radial  direction.  Let  the  flow  field  end  (at  radius  r  for  example),  and  let 
the  flow  properties  be  calculated  at  a  grid  node  at  this  radial  distance  r  and 
axial  distance  x+(l/2)Ax  (e.g.,  cross  in  Fig.  6).  The  previous  right  state  Sr 
(pr,  pr,  Mr,  £r)  at  the  coordinate  point  [x ,r-(l/2)Ar]  is  already  known.  Now, 
by  creating  a  pseudo  or  an  artificial  left  state  Sj  outside  the  flow  field  at 
the  coordinate  point  [x ,r+(l/2)Ar] ,  which  has  identical  flow  conditions  as  that 
of  the  right  state  (i.e.,  pj  =  pr,  =  pr ,  Mj  =  Mr,  and  =  £r) ,  the  Riemann 
problem  can  then  be  solved  in  the  usual  manner  to  obtain  the  required  flow  con¬ 
ditions  downstream  at  the  desired  coordinate  point  [x+(l/2)Ax ,r]  on  the  bound¬ 
ary.  Finally,  note  that  the  previous  specification  of  this  boundary  condition 
is  only  approximate,  and  minor  wave  reflections  can  sometimes  occur,  especially 
if  the  gradients  of  the  actual  flow  properties  across  the  boundary  are  relative¬ 
ly  large. 


Now  consider  the  case  of  a  rigid  surface  that  compresses  or  expands  the 
flow.  For  the  purpose  of  illustration,  let  the  flow  be  above  the  surface  of  a 
body.  In  this  case  the  left  state  Sj  (pj,  pj,  Mj,  ^1>  is  known  at  the  coordin¬ 
ate  point  [x,r+(l/2)Ar]  within  the  flow  (see  Fig.  6),  and  the  flow  properties 
for  the  downstream  state  at  coordinate  point  [x+(l/2) Ax ,r]  need  to  be  computed. 
Although  this  grid  node  may  be  located  inside  or  outside  of  the  flow,  its  state 
is  still  computed  for  either  case,  and  furthermore  employed  in  the  computations 
in  subsequent  steps.  In  order  to  obtain  the  flow  properties  at  this  downstream 
grid  node,  an  artificial  state  Sr  is  first  created  for  the  coordinate  point 
[x,r-(l/2)Ar]  that  is  located  inside  the  body. 

The  properties  pr,  pr,  Mr,  and  of  this  artificial  right  state  are 

obtained  by  applying  a  symmetric  boundary  condition.  An  artificial  flow  with  a 
right  wave  is  created  inside  the  body  to  simply  counteract  the  flow  outside  such 
that  the  flow  deflection  of  the  resulting  slip  stream  is  tangent  to  the  surface 
of  the  body.  From  the  point  of  view  of  symmetry,  therefore,  pr  =  p^,  pr  =  pj, 

M r  =  Mj,  and  £r  =  5l  ~  28,  where  6  is  the  flow  deflection  angle.  This  symmetry 
gives  left  and  right  waves  of  the  same  type  and  with  the  same  strength,  and  the 
angle  6  is  positive  for  a  left  shock  wave  and  negative  for  a  left  rarefaction 
wave.  The  deflection  angle  6  can  be  related  to  the  angle  of  the  body  surface  by 
8  =  -  n/2  +  tan-^(drjj/dx)  or  8  =  -^r  +  n/2  -  tan"^(dr^/dx)  ,  where  dr^/dx  is 

the  local  slope  t.an(6j,)  of  the  surface.  Hence,  £r  =  n  -  Si  -  2tan-1  (dr^/dx )  is 
the  final  result  for  the  flow  properties  at  the  artificial  state  Sr .  Now  that 
both  left  and  right  states  are  known,  the  Riemann  problem  can  be  solved  and 
random  sampling  applied  in  the  usual  manner  to  determine  the  new  flow  properties 
at  the  downstream  grid  node. 

Finally,  consider  the  case  of  a  free-jet  boundary.  For  illustration 
purposes  again,  let  the  free-jet  flow  be  above  or  to  the  left  of  the  boundary. 

As  in  the  previous  case,  the  left  state  Sj  (pj ,  pj,  Mj ,  is  known,  the  flow 

properties  at  the  downstream  grid  point  need  to  be  computed,  and  an  artificial 
right  state  Sr  needs  to  be  established.  A  symmetric  boundary  condition  is  again 
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applied,  with  a  slight  difference  in  the  last  step.  Firstly,  pr  =  pj,  pr  =  p1# 
Mr  =  M1#  and  £r  =  -  26,  where  6  is  the  again  flow  deflection  angle.  This 

symmetry  also  gives  left  and  right  waves  of  the  same  type  and  with  the  same 
strength,  and  the  angle  6  is  positive  for  a  left  shock  wave  and  negative  for  a 
left  rarefaction  wave.  The  flow  deflection  angle  can  again  be  expressed  in  the 
form  8  =  $1  ~  n/2  +  6j,  or  6  =  ~5r  +  n/2  -  6^ ,  where  8jj  is  the  angle  of  the  free- 
jet  surface.  However,  the  free-jet-surface  angle  6^  is  unknown  and  the  flow 
deflection  angle  6  cannot  be  calculated  from  6],.  This  angle  has  to  be  obtained 
from  known  boundary  information. 

The  information  required  is  that  the  pressure  at  the  free-jet  surface 
is  equal  to  the  pressure  of  the  fluid  or  gas  surrounding  the  free  jet.  Hence, 
the  strengths  p#/pi  and  p,/pc  of  the  left  and  right  waves  are  known,  and  the 
type  of  waves  are  also  known.  Based  on  the  known  strength  of  the  left  running 
wave  the  flow  deflection  angle  can  readily  be  calculated  for  a  left  shock  wave 
(p#/pi  >.1)  or  a  left  rarefaction  wave  (p*/pi  <  1),  by  using  the  appropriate 
equations  for  oblique  shock  or  rarefaction  waves.  In  the  present  case  this  is 
particularly  simple,  by  making  use  of  previously  introduced  equations.  The  flow 
deflection  angle  8  is  the  negative  of  the  previously  defined  function  X(n,M), 
that  is,  8  =  -  X(p#/pi,Mj)  in  our  case.  This  then  gives  the  flow  angle  at  the 
right  state  as  £r  =  kl  +  2X(p#/pi ,Mj ) .  This  completes  the  specification  of  the 
right  artificial  state,  so  that  the  normal  Riemann  problem  can  be  solved  and 
random  sampling  applied  to  get  all  the  flow  properties  at  the  downstream  state. 


2.8  Operator-Spl itt ina  Technique  to  Obtain  Axisvmmetric-Flow  Solutions 

The  previous  analysis  is  complete  for  the  random-choice  method  to  solve 
two-dimensional,  planar,  steady,  supersonic  flow  problems,  that  is  problems  for 
which  the  inhomogeneous  term  h(W,r)  in  Eqs .  1  and  2  is  nonexistent  or  just  zero 
(because  a  =  0) .  For  axisymmetric  flows,  however,  this  is  not  true  because  the 
inhomogeneous  term  is  nonzero  (a  =  1).  Some  additional  analysis  is  required  to 
correct  the  planar-flow  solution,  which  is  determined  first  for  the  case  of  axi- 
symme trie -flow  problems,  such  that  the  axisymmetric-f low  solution  of  interest  is 
then  obtained.  The  operator-splitting  technique,  originally  introduced  to  the 
random-choice  method  by  Sod^l  to  include  the  inhomogeneous  term  as  a  correction 
to  solve  one-dimensional  unsteady  flow  problems  with  an  area  change,  is  also 
used  in  this  study. 

In  the  case  of  the  operator-splitting  technique,  the  first  step  is  to 
solve  the  homogeneous  set  of  equations  in  conservation  form. 


|^tf (W)  ]  +  ^tgd)]  =  0, 


(42) 


for  W  =  [p,  pu,  pv,  e]  by  means  of  the  past  analysis  for  the  Riemann  problem. 
Once  this  Riemann  solution  V  for  a  planar  flow  is  known,  the  second  step  is  to 
use  the  set  of  equations  given  by 


|5-[f(W)]  =  h(W,r) 


(43) 


to  correct  this  planar  solution  W  to  obtain  the  axisymmetric-f low  solution  W. 
Because  the  right  hand  side  of  this  equation  is  known,  it  becomes  an  ordinary 


differential  equation,  and  the  correction  at  the  node  on  the  i**1  column  and  j1*1 
row  is  obtained  by  means  of  a  Cauchy-Euler  finite-difference  scheme  in  the  form 

f(Wj)  =  f(wj)  +  h(Wj.rj)  A*.  (44) 

where  Ax  is  the  incremental  distance  from  the  previous  left  and  right  states 
to  this  downstream  grid  node.  The  solution  W  can  then  be  written  in  the  form 


p  u 
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or  much  more  simply  as  pu  =  Xpu  =  Gj ,  p+pu2  =  p+Xpu2  =  G2 ,  puv  =  Xpuv  =  G3 ,  and 
u(e+p)  =  Xu(e+p)  =  G4,  where  X  =  1  -  (Ax/r)v/u  is  a  correction  factor.  If  X  is 
equal  to  unity  there  is  no  correction  required,  and  as  X  deviates  from  unity  the 
correction  becomes  less  accurate.  From  the  first  and  third  of  the  last  four  of 
these  equations  it  can  be  clearly  seen  that  the  radial  velocity  v  will  always 
be  invariant  during  the  correction  procedure. 


The  values  of  Gj ,  Gj ,  G3 ,  and  G4  and  can  be  manipulated  to  yield  values 
of  p,  p,  u,  v,  e,  M,  and  5 .  The  quadratic  equation 


with 


and 


Ap2  -  2Bp  -  C  =  0, 

A  =  [G3/G1  -  2G4]/Gi  =  v2-2(e+p)/p, 

B  =  -  [y/(y-1)]G2  =  -  iy! (y~D l(Xpu2  +  p) , 
C  =  [<y+1>/<y-1>1Gi  =  [(r+l)/(y-l)HXpu)2, 


(46) 


is  first  established  in  terms  of  the  density  p,  and  then  the  solution  for  the 
density  is  given  by 


p  =  [B  +  (B2  +  AC)1/2J/A.  (47) 

The  other  primitive  variables  then  follow  from  u  =  Gj/p  =  Xpu/p,  v  =  G3/G1  =  v, 
p  =  G2  -  G^/p  =  p  +  Xpu2(l-p/p),  and  e  =  PG4/G1  -  p  =  (e  +  p)p/p  -  p.  Finally, 
M  =  [(u2  +  v2)/(yp/p) ]2^2  and  £  =  n/2  -  tan-2(v/u),  completing  the  corrected 
axisymmetr ic-f low  solution  for  this  grid  node  (i,j). 

The  procedure  of  correcting  the  initial  planar-flow  solution  to  obtain 
the  axisymmetric-f low  solution  by  means  of  the  operator-splitting  technique  is 
explicit,  and  this  keeps  the  entire  random-choice  method  explicit.  Note  also 
that  the  boundary  conditions  are  not  used  in  the  operator-splitting  technique, 
because  they  enter  only  during  the  procedure  of  obtaining  the  initial  Riemann 
solution  for  the  planar  flow.  The  correction  scheme  by  means  of  the  operator¬ 
splitting  technique  is  applied  in  the  random-choice  method  at  each  grid  point 
within  the  flow  field,  except  for  nodes  on  the  origin  (r  =  0).  At  the  origin 
no  correction  is  needed,  because  the  correction  factor  X  =  1  -  (Ax/r)v/u  equals 
unity  by  virtue  of  the  fact  that  v/r  equals  zero  at  the  origin.  A  correction 
is  normally  done  just  after  the  initial  Riemann  problem  for  the  planar  flow  is 
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computed,  or  just  after  an  entire  column  is  computed.  Therefore,  each  initial 
Riemann  solution  for  the  planar  flow  is  obtained  from  both  the  left  and  right 
states  that  have  been  previously  corrected  to  give  the  axisymmetr ic-f low 
solution.  Making  the  correction  at  every  half-step  is  important  for  improving 
the  resolution  and  accuracy  of  the  numerical  results,  without  much  additional 
effort.  This  type  of  correction  scheme  is  first  order  and  can  be  called  an 
'asymmetric'  operating-splitting  technique^  ,5.*s(  an(j  is  use(j  most  often  in 
the  computations  in  this  report. 

There  is  also  another  first-order  equivalent  correction  scheme  called  a 
'symmetric'  operator-splitting  technique^, 58.  por  ^jj£s  scheme  the  equations 
are  the  same  as  those  already  given  in  this  section.  However,  in  this  scheme 
the  correction  is  applied  twice  in  succession  for  each  grid  point  in  each  alter¬ 
nate  column  (in  the  radial  direction)  and  not  at  those  columns  in  between.  This 
scheme  is  seldom  used  in  the  present  study,  although  it  gives  equivalent 
results  with  the  same  computational  effort. 

There  is  also  a  second-order  correction  technique  due  to  MacCormack®'^®. 
The  two  equations  for  this  predictor-corrector  technique  are 


f(wj) 


f(\?j) 


f<wj) 


+  h( Wj , r j )  Ax, 


(48) 


[  «*] 


)  +  f(W])  + 


h(Wj.rj)  Ax], 


where  W  is  the  first-order  intermediate  result  of  the  first  step  by  forward 
differencing  and  W  is  the  final  second-order  solution  from  the  second  step  by 
backward  differencing.  This  scheme  is  infrequently  used  in  this  study,  because 
it  require;,  more  computations  per  step  without  increasing  the  accuracy. 

A  second-order  correction  scheme  would  not  normally  be  used  with  a 
first-order  random-choice  method,  because  there  is  no  increase  in  accuracy  to 
be  expected.  However,  in  some  cases  numerical  noise  can  be  reduced,  and  this 
can  be  advantageous  and  justify  its  use.  This  is  not  the  case  in  the  present 
invest igat ion . 

It  is  worthwhile  mentioning  here  that  the  introduction  of  an  operator- 
splitting  technique  in  the  solution  of  axisymmetr ic-f low  problems  contributes 
numerical  error.  The  technique  gives  accurate  results  only  when  the  correction 
for  the  inhomogeneous  term  is  relatively  small.  In  axisymmetric  flows  for  which 
the  radius  r  includes  zero,  the  correction,  because  of  the  1/r  factor  in  h(W,r) , 
becomes  unduly  large  and  is  therefore  inaccurate.  Hence,  numerical  results  for 
some  parts  of  axisymmetric  flow  fields  which  include  small  values  of  radial  dis¬ 
tance  should  not  be  expected  to  be  computed  accurately.  Furthermore,  and  very 
importantly,  parts  of  axisymmetric  flow  fields  away  from  the  center  of  symmetry 
but  which  originate  or  depend  on  parts  of  the  flow  field>that  include  the  center 
of  symmetry  may  not  be  as  accurate  as  expected. 

Sod^  has  suggested  recently  that  more  accurate  numerical  results  can  be 
obtained  with  the  operator-splitting  technique  if  the  value  of  the  radius  r  in 
the  inhomogeneous  term  h(W,r)  is  computed  at  the  sampling  point,  rather  than  at 
the  location  of  the  downstream  grid  node.  Although  the  computational  effort  of 
adding  this  improvement  is  virtually  negligible,  most  RCM  codes  do  not  include 
this  improvement.  At  large  radii  or  numerous  grid  points  from  the  origin  the 
improvement  is,  of  course,  inconsequential.  For  small  radii  or  less  than  a  few 
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grid  nodes  from  the  origin  the  operator-splitting  correction  can  then  change. 
However .  these  changes  were  discovered  to  be  insignificant  in  the  solution  as  a 
whole,  and  have  not  been  included  in  most  of  the  numerical  results  to  be  given 
in  later  chapters. 


3.0  PRACTICAL  METHOD  OF  REDUCING  RANDOM-CHOICE-METHOD  COMPUTATIONS 

In  numerical  computations  of  planar  and  axisymmetric  steady  supersonic 
flows,  large  portions  of  the  flow  field  can  be  entirely  free  of  wave  motion  or 
have  insignificant  wave  motion.  These  portions  of  the  flow  field  need  not  be 
computed  in  most  cases,  which  can  result  in  a  marked  reduction  in  computational 
time  with  a  significant  saving  in  computational  cost. 

Consider  the  free-stream  flow  ahead  of  a  planar  or  axisymmetric  body,  as 
shown  in  Fig.  7.  Because  the  wave  emanating  from  the  body  is  swept  downstream 
as  it  moves  outward  in  the  supersonic  flow,  free-stream  conditions  ahead  of  the 
wave  exist  over  a  large  portion  of  the  flow  field  (labelled  A  in  the  figure). 
This  is  particularly  true  when  the  free-stream  Mach  number  is  large.  In  any 
case,  it  is  rather  easy  and  often  highly  beneficial  in  reducing  computational 
costs  to  write  the  computer  program  in  such  a  way  that  computer  computations 
are  skipped  at  all  of  these  grid  points  for  which  the  flow  properties  are  simply 
constant . 


Now  consider  the  flow  in  region  C  of  Fig.  7,  downstream  of  the  main  wave 
of  interest  in  region  B.  The  variation  in  the  flow  properties  throughout  region 
C  is  small  and  of  little  interest.  If  an  appropriate  boundary  condition  is  used 
along  the  dashed  line  shown,  the  flow  properties  in  region  C  can  be  eliminated. 
This  termination  line  must  be  sufficiently  far  behind  the  body  such  that  all  the 
important  or  interesting  parts  of  the  flow  field  are  included  in  the  numerical 
computations.  For  planar  and  axisymmetric  bodies  one  can  start  this  line  at  a 
radial  distance  equal  to  zero  and  an  axial  distance  of  approximately  one-half  to 
one  body  length  behind  the  end  of  the  body,  depending  on  the  flow  field.  Then, 
across  this  boundary  during  the  computations  the  flow  angle  can  be  taken  to  have 
a  linear  variation  (constant  gradient).  Hence,  the  flow  angle  at  the  artificial 
state  to  the  right  of  this  boundary  (r  is  given  by 

5j-l  =  2tj  -  tj+1  or  \t  =  -  5J+1,  (49) 

if  the  flow  is  above  the  boundary  and  the  flow  angle  £j+l  is  for  the  second  grid 
point  inside  the  terminal  line.  This  specification  of  £j-i  or  £r  for  the  right 
artificial  state  is  actually  equivalent  to  the  specification  of  a  solid  boundary 
with  a  surface  slope  n/2  -  2£j  +  £j+l  or  n/2  -  5j-l  somewhere  between  the  left 
and  right  states  that  straddle  this  terminal  line.  Therefore,  the  handling  of 
flow  computations  at  this  terminal  line  reduce?  to  the  simple  case  of  applying 
the  previous  solid-body  or  rigid-wall  boundary  condition  with  a  body  slope  given 
by  n/2  -  25 j  +  5j+1. 

As  the  flow  field  is  being  computed  from  the  terminal  boundary  line  to 
the  leading  shock  wave,  this  terminal  line  can  inserted  in  the  computations  by 
making  it  run  approximately  paralled  to  the  leading  shock  wave.  This  is  easily 
accomplished  by  selecting  a  convenient  number  of  grid  nodes  between  the  leading 
shock  wave  and  the  terminal  line.  Therefore,  the  numerical  computations  for  the 
main  outward  moving  wave  are  done  in  the  form  of  a  band  with  a  constant  width, 
which  is  labelled  B  in  Fig.  7. 
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4.0 


LIMITATIONS  OF  THE  RAM)OM-CHOICE  METHOD 


4 . 1  Steady  Supersonic  Planar  Flows 

The  random-choice  method  can  be  applied  to  obtain  solutions  to  numerous 
different  types  of  planar-flow  problems  (see  the  numerical  solutions  presented 
in  chapter  5) .  The  only  limitation  is  that  the  steady  flow  has  to  be  supersonic 
everywhere  in  the  computed  flow  field.  Hence,  for  flows  over  bodies,  oblique 
shock  waves  must  normally  be  attached  either  to  the  leading  edge  of  the  body  or 
one  of  its  protuberances.  Supersonic  flows  around  bodies  with  a  blunt  leading 
edge  or  supersonic  flows  that  contain  'pockets'  of  subsonic  flow  cannot,  there¬ 
fore,  be  computed. 

It  is  worthwhile  to  note  that  the  random-choice  method  cannot  calculate 
all  supersonic  flows  over  wedges  for  which  the  oblique  shock  is  attached  to  the 
leading  edge,  because  an  attached  shock  can  still  have  a  subsonic  flow  behind  it 
if  the  wedge  angle  is  sufficiently  large^®-^^.  Consequently,  the  limitation  in 
the  application  of  the  random-choice  method  for  solving  steady  planar  flows  is 
most  properly  that  the  flow  must  always  be  supersonic  and  not  the  criterion  of 
shock  attachment  or  detachment. 


4 .2  Steady  Supersonic  Axisymmetric  Flows 

The  random-choice  method  can  also  be  used  to  solve  many  different  types 
of  axisymmetr ic-f low  problems.  However,  its  application  is  more  limited  in  the 
case  of  axisymmetric  flows  than  planar  flows.  Besides  the  obvious  limitation 
that  the  steady  flow  has  to  be  supersonic  everywhere  within  the  computed  flow 
field,  there  is  one  additional  limitation  that  restricts  the  application  of  the 
random-choice  method,  and  also  one  additional  difficulty  that  can  make  the  grid 
spacing  undesirably  small. 

The  one  additional  limitation  stems  directly  from  using  a  first-order 
random-choice  method  for  which  the  planar  solution  is  computed  first  and  then 
corrected  by  means  of  the  operator-splitting  technique  to  obtain  the  desired 
axisymmetr ic-f low  solution.  For  certain  axisymmetric  flow  conditions  which 
involve  relatively  large  flow-deflection  or  turning  angles  through  shock  waves, 
the  corresponding  planar  solution  does  not  always  exist.  Hence,  there  is  no 
planar  solution  in  such  cases  to  correct  to  obtain  the  desired  axisymmetric-f low 
solution.  As  a  result  of  this  unfortunate  limitation,  axisymmetric  flows  with 
such  large  flow  deflections  cannot  be  computed,  or  flow  fields  cannot  be  fully 
computed . 

This  limitation  is  well  illustrated  by  a  plot  of  the  maximum  possible 
deflection  angle  through  an  oblique  shock  that  produces  a  sonic  flow  behind  it 
and  is  attached  to  a  wedge  and  a  cone.  Such  a  plot  as  shown  in  Fig.  8  can  be 
generated  from  the  tables  and  charts  in  the  NASA  Ames  hanlbook^^.  The  maximum 
deflection  angle  for  the  case  of  the  cone  is  substantially  larger  than  that  for 
the  wedge.  Consequently,  the  random-choice  method  cannot  solve  a  planar  Riemann 
problem  for  which  the  maximum  deflection  angle  lies  between  the  results  for  the 
cone  and  wedge,  because  the  oblique  shock  wave  for  the  case  of  the  wedge  would 
produce  a  subsonic  flow  downstream  of  the  shock.  For  these  conditions  the  shock 
would  normally  be  detached  from  the  wedge,  but  it  could  be  also  attached  with  a 
subsonic  flow,  as  mentioned  previously. 

The  previously  documented  limitation  can  be  rather  restrictive  for  the 
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computations  of  certain  types  of  axisymmetric-f low  problems.  For  example,  most 
military  projectiles  that  are  designed  for  flight  Mach  numbers  between  1  and  3 
have  noses  with  semivertex  angles  6  which  are  between  the  maximum  deflection 
angles  of  the  cone  and  wedge.  Furthermore,  they  generally  have  small  but  abrupt 
protuberances  in  their  body  shape  which  produce  large  flow  deflections.  Hence, 
the  flow  fields  around  such  bodies  cannot  be  computed  practically  by  the  present 
first-order  random-choice  method.  This  is  rather  unfortunate  because  a  lot  good 
experimental  data  exists  for  military  projectiles^*^-^,  which  would  have  been 
extremely  useful  in  verifying  the  capability  of  the  present  random-choice  method 
for  computing  axisymmetric-f low  fields. 

The  additional  difficulty  with  sometimes  having  to  use  an  undesirably 
fine  grid  mesh  when  computing  axisymmetric-f low  fields  (as  compared  to  planar 
ones)  is  now  described.  This  difficultly  arises  only  when  the  flow  field  is 
computed  at  the  origin  (i.e.,  r  =  0  or  r  =  0) ,  and  is  connected  directly  with 
the  use  of  a  first-order  random-choice  method  for  which  the  planar  solution  is 
computed  first  and  then  corrected  by  means  of  the  operator-splitting  technique 
to  obtain  the  axisymmetric-f low  solution.  In  this  case  the  correction  by  the 
operator-splitting  technique  for  a  small  radius  is  large  and  inaccurate,  as 
mentioned  at  the  end  of  the  previous  section. 

In  order  to  limit  this  type  of  inaccuracy  to  a  smaller  portion  of  the 
flow  field  near  the  origin,  a  finer  mesh  is  required,  because  this  type  of 
error  typically  diminishes  from  cell  to  cell  farther  from  the  origin  (and  not 
just  on  increasing  the  radial  distance)  .  The  inaccuracy  will  generally  become 
negligible  after  a  certain  number  of  cells  away  from  the  origin  (e.g.,  twenty). 
This  is  the  main  reason  why  a  finer  mesh  will  help  limit  the  spatial  extent  of 
the  error.  However,  a  finer  mesh  results  in  an  increase  in  computational  time 
and  cost. 

The  discussion  of  a  particular  example  can  helpful.  Consider  the  steady 
supersonic  flow  over  the  axisymmetric  body  shown  in  Fig.  7.  The  correction  by 
the  operator-splitting  technique  will  give  inaccurate  flow  properties  over  the 
entire  body  and  also  away  from  it  if  too  few  cells  are  used  (e.g.,  10  cells  per 
body  radius).  If  more  cells  are  used  the  flow  properties  will  be  inaccurate  in 
only  the  vicinity  of  the  nose  and  tail  of  the  body.  Because  of  the  inaccurately 
computed  flow  field  in  these  two  regions,  the  formation  of  both  the  front  and 
rear  shocks  and  their  subsequent  placement  in  the  flow  field  will  be  affected. 

If  this  type  of  error  extends  outward  for  about  20  cells,  then  100  cells  per 
body  radius,  say,  would  be  needed  to  try  limit  this  inaccurate  region  to  only  a 
small  part  of  the  flow  field.  The  use  of  100  cells  per  body  radius  is  generally 
undesirably  large  for  solving  such  problems,  as  compared  to  as  few  as  20  cells 
per  body  radius  for  a  wedge  flow  for  example.  This  example  helps  highlight  the 
unfortunate  problem  of  being  forced  to  resort  to  an  undesirably  fine  mesh  to 
retain  computational  accuracy  when  employing  the  present  random-choice  method 
for  solving  certain  axisymmetric-f low  problems  that  include  the  origin  in  the 
flow-field  solution. 


5.0  NUMERICAL  RESULTS  AND  DISCUSSION 


Numerical  results  from  the  random-choice  method  (RCM)  are  presented  and 
discussed  for  many  different  planar  and  axisymmetric  flow  problems,  in  order  to 
illustrate  the  capability  and  flexibility  of  the  RCM.  Results  are  given  first 
for  planar  flows  and  then  for  axisymmetric  flows.  In  each  case  the  solution  to 
the  simplest  problem  is  presented  first,  followed  by  problems  with  increasing 
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complexity.  In  some  cases 
numerical  or  RCM  results, 
be  that  a  perfect  diatomic 
7/5  or  1.40. 


theoretical  results  are  also  presented  to 
Furthermore,  for  each  problem  the  flow  is 
gas  or  perfect  air  having  a  specific-heat 
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5 . 1  Steady  Supersonic  Planar  Flows 

5.1.1  Flow  over  a  compressive  corner 

The  first  problem  involves  a  uniform  flow  at  a  flow  Mach  number  Mj  equal 
to  3  that  encounters  a  compressive  wedge  of  30  degrees.  This  wedge  deflects  the 
flow  by  means  of  a  single  oblique  shock  wave  attached  to  its  leading  edge.  The 
numerical  results  for  this  relatively  simple  problem  are  shown  in  Fig.  9,  where 
the  theoretical  oblique-shock-wave  location  is  also  depicted.  The  RCM  and  exact 
solutions  for  the  shock-wave  location  are  in  excellent  agreement,  although  the 
RCM  predicts  a  slightly  different  shock  angle.  On  the  average  the  random  posi¬ 
tion  of  the  shock  location  from  the  RCM  is  correct,  within  numerical  error.  The 
RCM  and  exact  solutions  for  the  uniform  flow  properties  behind  or  downstream  of 
the  oblique  shock  are  in  perfect  agreement,  as  should  be  expected  of  the  RCM 
for  this  simple  problem.  The  location  of  the  wedge  surface  in  the  RCM  computa¬ 
tions  is  also  shown  (for  interest).  Note  that  the  RCM  results  were  calculated 
with  a  grid  of  100  nodes  in  the  lateral  (y)  direction  and  100  in  the  longitudi¬ 
nal  (x)  direction.  However,  only  every  eighth  grid  node  in  the  longitudinal 
direction  is  shown,  and  this  number  is  indicated  in  the  inset  as  RCM  (8). 


5.1.2  Flow  over  an  expansive  corner 

The  second  problem  also  involves  an  initially  uniform  flow,  but  at  a 
flow  Mach  number  Mj  equal  to  2  over  an  expansive  corner  with  an  angle  of  30.53 
degrees.  The  flow  is  turned  by  a  centered  or  Prandtl-Meyer  rarefaction  wave 
that  emanates  from  the  sharp  corner,  as  depicted  in  Fig.  10a.  The  numerical 
(RCM)  and  exact  results  for  the  locations  of  the  head  and  tail  of  the  expansion 
fan  are  shown  in  the  figure  to  be  in  excellent  agreement.  Again,  the  RCM  and 
exact  results  for  the  flow  behind  the  rarefaction  wave  are  in  perfect  agree¬ 
ment,  as  should  be  expected  for  this  simple  problem.  Note  that  the  numerical 
results  have  been  computed  with  a  grid  with  100  nodes  in  the  y  direction  and 
100  nodes  in  the  x  direction,  but  results  are  given  for  only  every  eighth  node 
in  the  figure. 

RCM  and  exact  results  for  spatial  distributions  of  the  flow  properties 
through  the  rarefaction-wave  fan  are  shown  in  Figs.  10b  and  10c.  Profiles  in 
the  longitudinal  direction  x  (along  the  dashed  horizontal  line  in  Fig.  10a)  are 
given  first  in  Fig.  10b  and  profiles  in  the  lateral  direction  y  (dashed  vertical 
line  in  Fig.  10a)  appear  last.  In  each  case  the  agreement  between  the  RCM  and 
exact  results  is  excellent,  mainly  because  a  sufficient  number  of  nodes  was  used 
to  guarantee  good  agreement . 


5.1.3  Gradual  compressive  corner  that  gives  a  focussed  compression  wave 

The  problem  in  this  case  is  the  turning  of  an  initially  uniform  steady 
supersonic  flow  by  a  contoured  compressive  corner  that  produces  a  focussed 
compression  wave  which  results  in  an  oblique  shock  wave  (see  Figs.  11a  to  lid) 
This  type  of  problem  has  recently  been  studied  with  an  approximate  analysis  by 
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Emanue 1^2-63 #  The  focal  point  of  the  compression  wave  is  the  origin  of  the 
primary  oblique  shock  wave,  a  slip  stream  (dashed  line)  that  separates  the  flow 
through  the  shock  wave  from  the  isentropic  flow  through  the  oblique  compression 
wave,  and  a  weak  secondary  disturbance  that  can  be  either  an  oblique  shock  or  a 
Prandtl-Meyer  rarefaction  wave,  depending  on  the  initial  flow  conditions  and 
final  wall  angle.  The  secondary  wave  moves  outward  from  the  focal  point  and 
subsequently  reflects  first  from  the  wedge  surface,  then  interacts  with  the  slip 
stream  and  eventually  overtakes  the  primary  shock  wave.  The  free-stream  flow 
Mach  numbers  Mj  for  the  four  sets  of  numerical  results  presented  in  Figs.  11a 
to  lid  are  1.3,  1.64,  2.0,  to  4.0,  respectively,  as  indicated  in  the  figures. 
These  initial  conditions  correspond  exactly  to  some  of  those  used  in  Emanuel's 
investigations.  Note  that  the  flow  Mach  number  M  and  pressure  ratio  P  =  p/pj 
are  indicated  for  various  flow  regions  in  the  four  figures.  Mj  and  pj  are  the 
initial  or  free-stream  conditions. 

The  contoured  boundary  surface  to  produce  the  focussed  compression  wave 
is  not  a  known  function,  but  has  to  be  determined  for  the  computations^^.  In 
the  present  study  it  was  generated  numerically  with  the  RCM.  The  flow  is  first 
reversed  so  that  it  goes  through  a  centered  Prandtl-Meyer  expansion  wave  (e.g., 
from  M  =  1.178  to  Mj  =  1.64  in  Fig.  lib).  This  expansion  flow  is  computed  and 
these  computations  include  tracing  a  stream  line  at  some  reasonable  distance 
away  from  the  corner.  This  line  then  becomes  the  contoured  boundary  surface 
for  the  computations  of  the  reverse  flow  through  the  focussed  compression  wave. 

The  graphical  results  shown  in  Figs.  11a  to  lid  were  constructed  from 
the  numerical  or  RCM  results,  which  used  a  rather  fine  grid  of  200  nodes  in  the 
y  direction.  The  number  of  nodes  in  the  x  direction  varied,  being  574,  539,  496 
and  260  for  Figs.  11a  to  lid,  respectively.  None  of  the  results  in  the  figures 
are  theoretical,  although  exact  results  could  have  been  obtained  by  employing 
the  method  of  characteristics  and  Rank ine-Hugoniot  equations,  and  approximate 
predictions  are  available  also  from  the  work  of  Emanuel^^-63 #  Note  that  the 
fan  of  characteristics  for  the  focussed  compression  wave  has  been  constructed 
as  straight  lines,  as  would  occur  for  the  theoretical  solution  of  a  focussed 
compression  wave.  The  flow  properties  inside  the  fan  (such  as  the  pressure, 
density,  flow  Mach  number,  and  flow  angle)  were  then  checked  to  see  if  they 
were  constant  along  these  straight  lines.  They  were  indeed  very  close  to  being 
constant  along  these  lines,  showing  that  the  lines  are  indeed  characteristic 
lines.  It  was  also  found  that  the  focal  point  was  well  defined  in  a  very  small 
region  of  space  (shown  as  a  solid  dot),  corresponding  to  an  area  about  the  size 
of  4AxAy.  The  numerical  error  in  the  constancy  of  flow  properties  along  the 
characteristic  lines  and  in  defining  the  focal  point  is  small  because  of  the 
fine  grid  used  in  the  computations. 

Some  of  the  numerical  results,  such  as  the  oblique  shock  angle,  were 
compared  to  those  from  the  analyses  of  Emanue 1^2-63 f  an(j  the  agreement  was  very 
good,  provided  that  his  weak  shock  solutions  wr,re  used.  Such  comparisons  are 
not  presented  here,  mainly  because  his  analysis  is  not  exact  and  any  small 
differences  in  the  comparisons  are  therefore  not  known  to  result  from  the  RCM 
or  his  approximate  solution. 


5.1.4  Flow  over  a  symmetric  double-wedae 

Numerical  computations  of  a  supersonic  flow  (M^  =  2)  over  a  double  wedge 
are  presented  in  Fig.  12a  and  12b.  This  symmetric  double  wedge  is  1  unit  long, 
has  a  thickness  of  0.0525  units,  and  its  semivertex  angle  is  6  degrees.  The 
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number  of  grid  points  was  100  in  the  y  direction  and  220  in  the  x  direction. 

The  wave  diagram  appears  first  in  Fig.  12a.  It  shows  the  trajectories 
of  both  the  front  and  rear  shocks  from  the  leading  and  trailing  edges  of  double 
wedge,  as  well  as  showing  clearly  the  spreading  of  the  expansion  fan  emanating 
from  the  double-wedge  corner  formed  by  the  junction  of  the  bases  of  the  wedges. 
This  fan  eventually  expands  such  that  its  head  overtakes  the  front  shock  and  the 
tail  falls  into  the  rear  shock,  after  which  the  interaction  of  the  shocks  and 
expansion  wave  continually  causes  the  strengths  of  both  shocks  and  the  expansion 
wave  to  diminish.  The  flow  Mach  number  M  and  pressure  ratio  P  (p/pj)  in  regions 
having  a  constant  states  are  also  indicated. 

The  exact  solution  for  the  trajectories  of  the  front  and  rear  shocks  and 
the  head  and  tail  of  the  expansion  wave,  prior  to  their  interaction,  are  shown 

for  comparison.  The  agreement  in  general  is  very  good,  although  the  front  shock 

trajectory  from  the  ROM  lies  slightly  out  front  of  that  from  the  exact  solution. 

The  agreement  can  be  improved,  of  course,  by  using  more  grid  points. 

Pressure  distributions  through  the  entire  wave  produced  by  the  double 
wedge  are  shown  in  Fig.  12b  at  different  lateral  distances  away  from  the  double 
wedge.  These  profiles  show  clearly  how  the  wave  from  the  double  wedge  evolves 
into  a  decaying  N-shaped  wave  at  larger  lateral  distances. 


5.1.5  Flows  in  duct  inlets 


Supersonic  flows  into  ducts  of  changing  cross-section  always  involve 
oblique  shock  and  rarefaction  waves.  Numerical  solutions  for  seven  duct  inlet 
problems  are  now  given  and  discussed  (see  Figs.  13  to  19).  For  these  solutions 
the  number  of  grid  points  in  the  y  direction  was  100  for  the  first  three  cases 
and  150  for  the  remainder,  whereas  the  number  of  grid  nodes  in  the  x  direction 
was  350,  300,  500,  360,  290,  350,  and  350  for  the  seven  cases.  In  all  cases 
the  wave  diagrams  are  presented,  along  with  numerical  results  for  the  flow  Mach 
number  and  pressure  in  regions  of  constant  state.  Further,  exact  solutions  are 
given  for  the  trajectories  of  oblique  shock  waves,  rarefaction-wave  heads  and 
tails,  and  slip  streams.  In  the  seven  cases  the  RCM  and  exact  solutions  are  in 
good  agreement,  and  this  agreement  can  be  improved  by  using  more  grid  nodes. 

The  first  case  is  for  a  uniform  free-stream  flow  (Mj  =  2.0)  entering  a 
duct  inlet  for  which  the  bottom  duct  surface  is  a  compressive  wedge  of  6  degrees 
and  the  upper  duct  surface  is  a  flat  plate  parallel  to  the  free-stream  flow  (see 
Fig.  13).  The  oblique  shock  wave  emanating  from  the  lower  wedge  crisscrosses 
the  duct  four  times  as  it  reflects  alternately  from  the  upper  and  lower  duct 
walls.  Its  last  reflection  from  the  bottom  surface  is  a  Mach  wave,  because  at 
this  last  reflection  point  the  wedge  surface  is  simply  terminated  by  a  corner 
to  make  the  duct  area  constant  thereafter. 

♦ 

The  second  case  is  similar  to  the  first  one  but  has  two  distinct  dif¬ 
ferences  in  the  inlet-duct  geometry.  The  bottom  wedge  surface  with  an  angle  of 
6  degrees  is  no  longer  terminated  by  a  corner,  and  the  upper  surface  no  longer 
is  a  flat  plate  but  instead  is  a  compressive  wedge  of  2  degrees  (see  Fig.  14). 
Computations  of  this  flow  field  has  been  terminated  just  ahead  of  the  second 
oblique  shock-wave  reflection  at  the  upper  surface,  because  the  next  reflected 
shock  has  subsonic  flow  behind  it.  Hence,  any  additional  flow  field  cannot  be 
solved  by  means  of  the  random-choice  method,  illustrating  one  of  its  primary 
limitations  that  was  mentioned  previously. 
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For  the  third  case  of  an  inlet-duct  flow  the  upper  duct  surface  is  a 
flat  plate  aligned  with  the  free-stream  flow  (Mj  =  2.0)  and  the  lower  duct  sur¬ 
face  is  also  a  flat  plate  but  with  an  expansive  corner  at  the  leading  edge  (see 
Fig.  IS).  In  this  case  a  Prandtl-Meyer  expansion  wave  emanates  from  this  sharp 
corner  and  continues  to  spread  as  it  crisscrosses  the  duct  through  reflections 
from  the  upper  and  lower  surfaces.  Note  that  in  this  case  the  flow  is  expanded 
to  a  lower  pressure  and  higher  Mach  number  by  the  crisscrossing  rarefaction  wave 
such  that  it  adjusts  to  the  larger  flow  area  in  the  duct.  In  the  first  two 
cases  the  flow  was  compressed  to  a  higher  pressure  and  lower  Mach  number  by  the 
crisscrossing  shock  wave  in  the  process  of  adjusting  the  flow  to  the  smaller 
duct  area  . 

In  the  fourth  case  the  upper  and  lower  surfaces  of  the  inlet  duct  are 
both  compressive  wedges  (Fig.  16).  The  upper  wedge  has  an  angle  of  6  degrees 
and  the  lower  one  has  an  angle  of  only  3  degrees.  In  this  case  oblique  shock 
waves  originate  from  the  leading  edges  of  the  wedges  and  crisscross  the  duct. 
Their  interactions  produce  slip  streams,  and  the  subsequent  interaction  of  the 
oblique  shock  waves  and  slip  streams  produce  additional  weak  reflected  waves. 

In  all  cases  the  flow  properties  in  various  regions  of  constant  states  can  be 
determined  accurately  to  machine  precision,  although  the  boundaries  of  these 
regions  are  less  accurate,  within  numerical  error. 

The  fifth  case  is  similar  to  the  last  one  with  two  compressive  wedges. 

In  this  case,  however,  the  compressive  wedges  are  both  6  degrees  and  the  two 
oblique  shock  waves  that  originate  at  their  leading  edges  therefore  have  the 
same  strength  (Fig.  17).  Consequently,  their  subsequent  interactions  before 
and  after  reflecting  at  the  side  walls  produce  slip  streams  with  identical  flow 
properties  on  each  side.  Hence,  slip  streams  do  not  appear  in  the  figure. 

In  the  sixth  case  the  upper  surface  of  the  inlet  duct  is  an  expansive 
wedge  of  2  degrees  and  the  lower  surface  is  a  compressive  wedge  of  6  degrees 
(see  Fig.  18).  The  main  features  of  this  flow  are  the  interaction  of  the  two 
oblique  shock  and  rarefaction  waves  which  produce  a  slip  stream  of  finite  width, 
that  is  a  slip  region.  Such  regions  are  routinely  handled  by  the  RCM. 

The  seventh  and  final  case  is  for  the  supersonic  flow  in  a  duct  inlet 
consisting  of  two  expansive  wedges,  one  of  1  degree  for  the  upper  surface  and 
the  other  of  2  degrees  for  the  lower  surface  (Fig.  19).  In  this  case  the  main 
features  are  two  Prandtl-Meyer  expansion  waves  that  emanate  from  the  two  sharp 
expansive  corners  and  crisscross  the  duct.  A  slip  stream  or  region  from  the 
interaction  of  the  two  rarefaction  waves  does  not  appear  because  it  does  not 
exist  for  such  an  isentropic  flow. 


5.1.6  Flow  over  a  parabolic  shaped  airfoil 

Numerical  computations  were  done  for  a  supersonic  flow  (Mj  =  2.0)  over 
a  parabolic  shaped  airfoil  with  an  arbitrary  length  L  and  maximum  thickness  to 
length  ratio  t/L  of  0.10.  The  symmetric  airfoil  shape  is  given  by  the  equation 
yj,/L  =  2(t/L)[x/L  -  (x/L)^].  For  this  particular  case  the  trajectories  of  the 
front  and  rear  shock  waves  look  like  those  sketched  in  Fig.  7,  and  overpressure 
signatures  as  a  function  of  longitudinal  distance  (x)  are  presented  in  Fig.  20. 
The  first  signature  is  for  the  overpressure  (p-pj)/ Pi  along  the  airfoil  surface, 
and  the  second,  third,  and  fourth  overpressure  signatures  are  at  lateral  distan¬ 
ces  (y)  from  the  airfoil  of  one,  ten,  and  one  hundred  body  thicknesses.  (These 
distances  are  indicated  in  the  figure  by  y  =  It,  lOt ,  and  lOOt . ) 
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The  evolution  of  the  wave  with  lateral  distance  into  an  N  shape  is 
readily  apparent.  One  can  also  observe  that  the  wave  is  decreasing  nonlinearly 
in  amplitude  and  simultaneously  stretching.  The  peak  overpressure  decays  almost 
as  y~l/2  at  larger  lateral  distances,  in  direct  areement  with  nonlinear  acoustic 
or  aerodynamic  theory  for  weak  but  finite-amplitude  waves  (with  nonlinear  wave 
steepening  effects) .  *  Note  that  linear  aerodynamic  theory  would  predict  that 
the  wave  amplitude  and  length  would  simply  be  constant  with  increasing  lateral 
d is  tance . 


The  computations  for  the  parabolic  shaped  airfoil  were  done  initially 
for  the  flow  around  the  airfoil  (0  <  y  <  1.9)  with  190  grid  nodes.  Then  the 
computations  were  continued  with  a  coarser  mesh  containing  only  every  second 
point,  leaving  only  about  one  hundred  grid  nodes  in  the  vertical  direction 
through  the  outward  moving  wave,  as  described  in  chapter  3.  The  first  three 
overpressure  signatures  were  computed  in  this  manner,  and  the  number  of  nodes 
was  sufficient  to  ensure  that  these  signatures  are  of  high  quality.  For  the 
fourth  one,  the  computations  were  redone  with  a  much  coarser  mesh,  that  is,  175 
grid  nodes  in  the  range  0  <  y  <  3.5.  Then  the  computations  were  continued  with 
an  even  coarser  mesh  containing  only  every  fifth  point,  thereby  leaving  only 
about  forty  grid  nodes  in  the  vertical  direction.  This  is  why  the  quality  of 
the  last  signature  is  not  quite  as  good  as  the  previous  three. 


5.1.7  Free-jet  flows  from  a  slot 

Two  sets  of  RCM  computations  have  been  completed  for  supersonic  streams 
that  issue  from  a  slot  into  a  surrounding  quiescent  fluid  of  constant  pressure 
(see  Figs.  21  and  22).  In  both  cases  the  flow  prior  to  emerging  from  the  slot 
with  parallel  walls  was  uniform.  The  resulting  free-jet  flows  differ  markedly 
in  structure  mainly  because  the  first  free  jet  (Fig.  21)  has  an  initially  higher 
pressure  than  that  of  the  surrounding  fluid  and  the  other  free  jet  (Fig.  22)  has 
a  lower  pressure.  In  both  cases  only  the  lower  half  of  the  free  jet  has  been 
computed  by  the  RCM;  the  upper  symmetric  half  has  been  added  as  dashed  lines  for 
completeness.  The  flow  Mach  number  M  and  pressure  ratio  P  (p/py)  are  also  shown 
for  various  regions  of  constant  state. 

In  the  first  case  of  a  free-jet  flow  from  the  slot  the  initial  Mach  num¬ 
ber  Mi  is  1.40  and  the  ratio  of  the  quiescent  fluid  pressure  to  that  initially 
in  the  flow  in  the  slot  is  0.749  (Fig.  21).  In  these  numerical  results  it  can 
be  readily  seen  that  symmetric  rarefaction  waves  emanate  from  the  slot  corners 
and  expand  the  flow  before  they  eventually  reflect  at  the  jet  boundaries  and 
form  compression  waves.  These  compression  waves  then  compress  the  free  jet  as 
they  come  to  a  focus,  eventually  back  to  the  original  flow  conditions  and  cross- 
sectional  area  at  the  slot  (within  numerical  error).  In  the  case  of  perfectly 
accurate  numerical  computations  (no  numerical  error),  this  process  would  be 
repeated  indefinitely  for  subsequent  free-jet  lobes. 

One  measure  of  the  accuracy  of  the  RCM  for  this  first  free-jet  computa¬ 
tion  is  the  degree  of  focussing  of  the  compression  waves  at  the  end  of  the  first 
lobe.  Although  the  RCM  does  not  produce  a  perfect  focal  point,  the  focus  was 
found  to  be  within  an  area  equal  to  4AxAy,  which  is  highly  commendable  of  the 
RCM  method.  Another  measure  of  the  accuracy  is  the  straightness  of  the  head  and 
tails  of  the  computed  rarefaction  and  compression  waves.  These  are  indeed  quite 
straight  because  sufficient  grid  points  were  used  in  the  computations  to  ensure 
this.  For  example,  the  number  of  grid  points  was  70  for  half  of  the  free  jet 
in  the  y  direction  and  350  in  the  x  direction. 
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For  the  second  case  of  a  free-jet  flow  the  initial  Mach  number  is 
1.60  and  the  ratio  of  the  quiescent  fluid  pressure  to  that  initially  in  the 
flow  in  the  slot  is  1.344  (Fig.  22).  In  these  numerical  results  it  can  be  seen 
readily  that  symmetric  oblique  shock  waves  emanate  from  the  slot  corners  and 
compress  the  jet  flow  before  they  eventually  reflect  at  the  jet  boundaries  as 
Prandt 1-Meyer  expansion  fans.  From  this  reflection  point  onward,  the  free  jet 
is  of  the  same  structure  as  for  the  previous  case.  Successive  lobes  containing 
alternate  expansion  and  compression  waves  are  simply  repeated,  as  if  a  slot  was 
located  at  the  end  of  the  shock-wave  portion  with  this  contracted  jet  area  and 
initial  jet  flow  properties  of  M  =  1.166  and  P  =  1.822.  As  for  the  last  case, 
the  RCM  is  very  successful  in  obtaining  these  numerical  results.  Note  that  the 
number  of  grid  points  in  the  y  direction  is  50  for  the  one-half  of  the  free  jet 
and  425  in  the  x  direction. 


5 .2  Steady  Supersonic  Axisvmmetric  Flows 
5.2.1  Flow  over  a  cone 


The  first  problem  involves  a  uniform  flow  at  a  flow  Mach  number  equal 
to  2.5787  that  encounters  an  infinitely  long  cone  with  a  semivertex  angle  0C  of 
15  degrees.  The  cone  axis  and  flow  are  aligned  (Fig.  23).  The  flow  is  partly 
turned  suddenly  through  an  oblique  shock  and  then  turned  gradually  thereafter, 
such  that  the  flow  eventually  becomes  tangent  to  the  cone  surface.  Furthermore, 
from  theoretical  considerations  it  is  well  known  that  the  shock  angle  os  is  con¬ 
stant  and  the  flow  properties  along  radial  lines  with  angle  ©  between  the  shock 
and  cone  surface  are  also  constant  *14-16 ,60-61 

The  flow  field  was  computed  with  a  square  mesh  (Ar  =  Ax)  and  200  grid 
nodes  in  the  radial  ( r)  direction.  RCM  computed  flow  properties  for  the  density 
ratio  p/pj,  pressure  ratio  p/pj,  flow  Mach  number  M,  and  flow  angle  £  are  shown 
in  Figs.  24a,  24b,  and  24c,  where  they  are  compared  to  exact  tabular  solutions^ 
(actually  very  accurate  numerical  solutions  of  ordinary  differential  equations). 
These  results  are  plotted  as  a  function  of  the  angle  0  from  the  axis  of  the  cone 
to  a  radial  line  in  the  flow  field,  at  three  different  distances  of  80,  160,  and 
320  grid  nodes  from  the  cone  apex.  Note  that  the  cone  surface  lies  at  an  angle 
of  15  degrees  and  the  exact  shock  location  is  27.85  degrees. 

The  agreement  between  the  RCM  and  exact  results  is  very  good  in  general. 
The  agreement  becomes  better  for  results  farther  downstream,  that  is  after  more 
grid  nodes  in  the  axial  direction.  One  reason  for  this  is  that  the  number  of 
nodes  in  the  radial  direction  through  the  RCM  computed  flow  field  increases  with 
distance  from  the  cone  apex.  The  data  in  Fig.  24b  has  twice  as  many  grid  nodes 
as  those  given  in  Fig.  24a,  and  the  data  in  Fig.  24c  has  twice  as  many  as  those 
in  given  Fig.  24b.  A  second  reason  is  that  the  flow-field  data  computed  farther 
away  from  the  cone  axis  is  more  accurate  by  the  RCM,  because  the  correction  of 
the  planar-flow  solution  to  obtain  the  axi-ymmetric-f low  solution  becomes  smal¬ 
ler  and  therefore  more  accurate. 

The  ax  is ymme trie- flow  solution  in  the  vicinity  of  the  cone  apex  can  be 
highly  inaccurate,  because  the  operator-splitting  correction  near  the  cone  axis 
is  large  and  incorrect.  The  pressure  ratio  across  the  shock  can  be  in  error  by 
as  much  as  70%.  However,  such  perturbations  were  found  to  decay  with  radial 
distance  to  almost  negligible  size  by  about  25  grid  zones  away  from  the  axis  of 
the  cone.  This  gives  a  good  indication  of  the  size  of  the  flow  field  affected, 
and  also  how  many  grid  points  are  required  in  this  and  other  ax isymmet ric-f low 
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problems  involving  cone-tipped  and  pointed  bodies,  to  make  sure  that  this  inac¬ 
curate  region  is  of  negligible  size  compared  to  the  entire  flow  field. 

A  close  inspection  of  the  RCM  and  exact  results  in  Figs.  24a  to  24c  will 
show  that  the  oblique  shock  location  from  the  RCM  is  out  front  (larger  9)  of  the 
exact  location.  The  difference  is  largest  for  the  data  given  in  Fig.  24a,  or 
nearer  the  cone  apex,  and  becomes  smaller  for  data  given  in  Figs.  24b  and  24c, 
or  as  the  distance  or  number  of  grid  nodes  from  the  cone  axis  increases.  The 
reasons  for  this  difference  in  shock  location  are  two-fold.  Firstly,  because 
the  operator-splitting  correction  can  result  in  flow-field  pressures  and  oblique 
shock-wave  strengths  that  are  too  high  in  the  vicinity  of  the  cone  apex,  the  RCM 
predicted  propagation  velocity  of  the  oblique  shock  wave  is  too  high  near  the 
cone  axis.  Secondly,  the  propagation  of  the  shock  depends  to  a  large  extent  on 
the  random-sampling  procedure.  This  procedure,  which  assigns  some  state  to  the 
next  grid  point  downstream,  is  always  completed  on  the  intermediate  planar  solu¬ 
tion,  and  then  this  downstream  state  is  corrected  to  give  the  final  solution  for 
axisymmetric  flow.  Because  the  oblique  shock  is  always  stronger  for  a  planar 
flow  than  an  axisymmetric  one,  the  sampling  procedure  on  the  intermediate  planar 
flow  solution  therefore  propagates  the  oblique  shock  wave  outwards  too  quickly, 
especially  in  the  vicinity  of  the  cone  apex.  This  effect  decays  with  increasing 
radial  distance,  but,  unfortunately,  it  usually  takes  40  to  50  grid  nodes  in  the 
radial  direction  before  this  effect  disappears. 

The  computation  of  an  oblique  shock  flow  fields  over  a  cone  is  an  excel¬ 
lent  test  of  the  capability  of  the  RCM  for  predicting  axisymmet ic-f low  fields, 
and  also  for  testing  out  an  initial  computer  program. 


5.2.2  Flow  over  a  cone-cvlinder 


Numerical  computations  were  done  for  the  classical  case  of  the  flow  over 
an  axisymmetric  body  with  a  conical  nose  on  an  infinite-length  cylinder.  The 
semivertex  angle  of  the  conical  nose  was  11.5  degrees.  From  these  computations 
the  pressure  coefficient  Cp  =  2(p  -  on  the  surface  of  the  body  was 

obtained,  where  p  is  the  static  pressure  on  the  body  surface.  These  RCM  results 
are  compared  to  an  exact  calculation^^  (method  of  characteristics)  in  Fig.  25. 


The  RCM  and  exact  results  are  in  good  agreement,  with  some  scatter  in 
the  RCM  results  near  the  cone  apex.  This  scatter  or  inaccuracy  is  due  mainly 
to  the  inaccuracy  of  the  operator-splitting  correction  near  the  cone  axis,  both 
of  which  were  just  discussed  for  the  previous  example.  Note  that  if  fewer  grid 
points  than  1100  in  the  x  direction  were  utilized  in  the  RCM  solution,  the  RCM 
results  would  then  become  worse  near  the  cone  apex. 


5.2.3  Flow  over  a  parabolic-nosed  cylinder 


Numerical  results  were  obtained  for  the  flow  over  a  cylindrical  body 
with  a  parabolic  nose,  for  three  different  free-stream  Mach  numbers  of  2.0,  3.0, 
and  3.92.  The  radius  R  of  the  infinite  length  cylinder  is  0.5  arbitrary  units, 
and  the  equation  for  its  parabolic  nose  having  a  length  L  of  3.8  units  is  given 
by  r^/L  =  2(R/L)[(x/L)  -  (x/L)^] .  Once  again,  the  pressure  coefficient  from  the 
RCM  data  is  compared  to  the  exact  results  (see  Figs.  26a,  26b,  and  26c),  and  the 
exact  results  are  from  the  method  of  characteristics^^.  For  the  three  different 
free-stream  Mach  numbers  the  agreement  between  the  RCM  and  exact  results  is,  in 
general,  good,  just  as  in  the  previous  example.  As  in  the  previous  case,  the 
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RCM  results  are  less  accurate  nearer  to  the  apex  of  the  parabolic  nose,  for  the 
same  reasons  given  previously.  In  each  case  the  number  of  grid  nodes  in  the  x 
direction  was  840,  and  the  mesh  used  was  square  (i.e..  Ax  =  Ar) . 

5.2.4  Flow  over  a  symmetric  double-cone 

The  RCM  was  used  to  compute  the  flow  field  for  the  supersonic  flow  over 
a  symmetric  double-cone.  The  cones  which  were  joined  at  their  common  bases  had 
a  semivertex  angle  of  10  degrees,  and  the  free-stream  Mach  number  M^  was  2.075. 
For  the  computations  a  square  mesh  was  employed,  having  100  grid  nodes  per  unit 
length  of  the  double-cone. 

The  RCM  solution  in  terms  of  axial  distributions  of  pressure  (P  =  p/pj) 
at  different  radial  distances  from  the  body  are  shown  in  Fig.  27.  These  results 
show  clearly  how  the  double-cone  produces  an  outwards  moving  wave  that  changes 
in  shape  and  eventually  evolves  into  an  N-shaped  wave.  Furthermore,  the  numeri¬ 
cal  results  also  illustrate  the  capability  of  the  RCM  for  obtaining  solutions  to 
such  problems . 

A  few  comments  are  worth  making  regarding  the  numerical  solution.  In 
the  first  few  pressure  distributions  near  the  double-cone  one  can  see  clearly 
the  initial  pressure  ramp  just  behind  the  front  shock  (which  originates  from 
the  conical  solution).  It  is  continually  being  eroded  by  the  rarefaction  wave 
that  emanates  from  expansive  corner,  as  it  overtakes  the  front  shock,  afterwhich 
the  rarefaction  wave  interacts  with  the  front  shock  wave  and  continually  reduces 
its  strength  with  increasing  radial  distance. 

The  flow  right  next  to  the  surface  of  the  double-cone  must  follow  the 
surface.  When  this  flow  turns  suddenly  at  the  nose  and  tail  of  the  body  the 
front  and  rear  shocks  are  formed.  A  little  farther  away  from  the  body  and  just 
after  the  front  shock,  the  flow  on  turning  suddenly  through  the  front  shock 
still  needs  to  be  turned  gradually  but  further  to  align  itself  with  the  body 
surface.  This  is  the  reason  for  the  pressure  ramp  following  the  front  shock. 

The  rarefaction  wave  signals  the  presence  of  the  corner  and  turns  the  flow  in  an 
attempt  to  realign  it  with  the  surface  of  the  second  cone.  In  the  process  it 
overturns  the  flow  toward  the  surface,  and  the  flow  then  has  to  gradually  turn 
away  from  the  surface.  This  produces  the  compression  wave  that  is  clearly  seen 
in  the  first  pressure  profile  just  ahead  of  the  rear  shock.  With  increasing 
radial  distance  from  the  body,  this  compression  wave  steepens  and  is  overtaken 
by  the  rear  shock.  This  combining  of  the  two  waves  is  the  reason  why  the  rear 
shock  wave  is  temporarily  stronger  than  the  front  shock.  Farther  away  from  the 
body  the  two  shock  waves  become  equal  in  strength  and  the  change  in  pressure 
between  them  becomes  linear. 

In  the  computations  for  the  results  given  in  Fig.  27,  the  free-stream 
flow  field  ahead  of  the  front  shock  was  not  computed  (see  chapter  3).  This 
reduces  the  computational  time  and  cost  by  about  10%.  However,  the  entire  flow 
field  behind  the  rear  shock  was  computed,  in  spite  of  the  fact  that  it  is  of 
little  interest.  Most  of  it  does  not  need  to  be  computed  and  can  be  eliminated 
as  described  in  chapter  3  to  reduce  computational  time  and  cost.  The  results  of 
additional  computations  for  which  this  lower  right  part  of  the  flow  field  was 
not  included  are  presented  in  Fig.  28.  The  computed  results  are,  of  course,  in 
perfect  agreement  with  those  in  Fig.  27,  and  the  extra  saving  in  computational 
time  and  cost  is  about  60%.  This  illustrates  clearly  the  benefit  to  be  gained 
by  eliminating  unnecessary  computations  of  certain  parts  of  the  flow  field. 
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The  first  pressure  profile  nearest  the  double-cone  in  Figs.  27  and  28 
contain  some  jaggedness.  This  is  due  to  inaccuracies  in  the  operator-splitting 
correction  near  the  double-cone  axis  for  axisymmetric-f low  computations.  These 
errors  can  be  reduced,  of  course,  by  increasing  the  number  of  grid  nodes  in  the 
RCM  computations,  especially  for  the  initial  flow  field  near  the  double-cone. 

To  illustrate  this  improvement,  the  initial  flow  field  in  the  region  contained 
by  the  lines  r  =  1.3  and  x  =  2.0  was  computed  first  with  twice  as  many  grid 
nodes  (200  nodes  per  unit  body  length).  Then  the  flow  field  containing  just  the 
outward  moving  wave  was  continued  outward  with  only  half  this  number  (100  nodes 
per  unit  body  length),  because  more  grid  nodes  are  not  needed  farther  away  from 
the  body.  The  new  results  are  shown  in  Fig.  29.  The  improvement  in  these  RCM 
results  over  the  previous  ones  in  Fig.  28  is  quite  obvious,  but  the  extra  compu¬ 
tational  time  and  cost  were  40%. 


5.2.5  Flow  over  a  parablolic-spindle  shaped  body 

Numerical  computations  were  done  for  a  supersonic  flow  (M^  =  2.0)  over 
a  parabolic-spindle  shaped  body  with  an  arbitrary  length  L  and  maximum  diameter 
to  length  ratio  D/L  of  0.10.  The  shape  of  the  axisymmetric  body  is  given  by  the 
equation  r^/L  =  2(D/L)[x/L  -  (x/L)^].  For  this  example  the  trajectories  of  the 
front  and  rear  shock  waves  look  like  those  sketched  in  Fig.  7,  and  overpressure 
signatures  as  a  function  of  axial  distance  (x)  are  presented  in  Fig.  30.  The 
three  signatures  are  for  the  overpressure  (p~pi)/pi  at  radial  distances  (r)  from 
the  body  of  one,  ten,  and  one  hundred  body  diameters. 

These  numerical  results  show  how  the  wave  evolves  with  radial  distance 
and  decays  into  an  N  shape.  The  first  two  signatures  are  in  good  agreement  with 
the  prediction  of  approximate  but  accurate  theory  by  Lighthill  and  Whitham  (see 
Ref.  24).  The  amplitude  of  the  third  signature  is  too  high  by  a  factor  of  about 
two  and  the  negative  overpressure  phase  is  too  large.  These  are  simply  due  to 
round-off  errors  from  multiple  single-precision  computations.  For  example,  the 
radial  flow  velocity  v  at  large  radii  is  very  small  such  that  the  flow  angle  £ 
measured  from  the  radial  axis  loses  significant  digits  and  becomes  virtually 
n/2.  This  particular  problem  could  be  removed  easily  by  simply  redefining  the 
flow  angle  to  be  measured  more  physically  from  the  axial  direction.  However, 
this  problem  does  illustrate  the  important  point  that  numerical  results  can 
become  saturated  with  roundoff  errors  in  the  case  of  weak  waves. 

Note  that  the  computations  for  the  parabolic-spindle  body  were  done 
initially  for  the  flow  around  the  body  (0  <  r  <  0.92  and  0  <  x  <  1.5)  with  920 
grid  nodes  in  the  radial  direction  and  1500  in  the  axial  direction.  Then  the 
computations  were  continued  with  a  coarser  mesh  containing  only  every  fifth 
point,  leaving  only  185  grid  nodes  in  the  radial  direction  through  the  outward 
moving  wave,  as  described  in  chapter  3. 


5.2.6  Free-iet  flows  from  a  circular  orifice 

The  first  numerical  computations  of  an  axisymmetric  free-jet  flow  are 
for  the  case  of  an  exit  flow  Mach  number  of  1.2  and  a  pressure  ratio  P  (p/pj) 
of  unity.  Because  the  pressure  of  the  fluid  outside  the  free  jet  is  only  one 
tenth  of  that  in  the  jet  at  the  exit,  the  free  jet  expands  on  emerging  from  the 
orifice.  Numerical  results  showing  this  expansion  by  tracing  the  jet  boundary 
are  presented  in  Fig.  31.  Also  shown  is  the  trajectory  of  the  oblique  shock 
wave,  commonly  called  the  barrel  shock,  which  starts  as  a  coalescence  of  Mach 
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wave  just  downstream  of  the  orifice  and  grows  in  strength  up  to  the  Mach  disc. 
The  Mach  disc,  and  the  Mach  reflection  near  the  periphery  of  the  Mach  disc,  are 
sketched  in  the  diagram,  because  they  cannot  be  computed  by  the  RCM  because  the 
flows  behind  them  are  subsonic. 

The  properties  of  the  flow  inside  the  free  jet  can  be  studied  easily 
with  the  RCM.  In  order  to  illustrate  some  of  the  features,  pressure  and  Mach 
number  distributions  with  radial  distance  are  included  as  insets,  one  in  the  top 
diagram  and  the  other  in  the  bottom  one.  The  changes  in  pressure  and  flow  Mach 
number  from  inside  the  barrel  shock  to  outside  the  free-jet  boundary  are,  for 
example,  rather  pronounced. 

The  second  numerical  computations  of  an  axisymmetric  free-jet  flow  are 
for  the  case  of  an  exit  flow  Mach  number  Mj  of  3.0  and  a  pressure  ratio  P  (p/pj) 
of  unity.  Because  the  pressure  of  the  fluid  outside  the  free  jet  is  now  twice 
as  high  as  that  in  the  jet  at  the  exit,  the  free  jet  contracts  on  emerging  from 
the  orifice,  as  shown  in  Fig.  32.  Besides  the  contracting  free-jet  boundary  the 
trajectory  of  the  inward  moving  oblique  shock  wave  is  presented.  Inside  this 
oblique  shock  the  flow  properties  are  simply  the  exit  conditions.  Across  this 
shock  the  flow  properties  change  drastically,  as  shown  by  the  insets  of  pressure 
and  flow  Mach  number,  and  then  more  slowly  to  the  jet  boundary. 

The  inward  moving  oblique  shock  eventually  nears  the  axis  of  symmetry, 
where  it  is  reflected.  Depending  on  the  flow  conditions  the  reflection  can  be 
either  a  regular  or  Mach  reflection.  In  the  present  case  it  appears  to  be  a 
regular  reflection,  although  the  Mach  stem  could  exist  and  be  rather  small.  The 
reflected  wave  is  simply  sketched  in  the  diagram,  because  the  flow  becomes  sub¬ 
sonic  and  the  computations  could  not  be  continued. 

Note  that  in  solving  the  two  free-jet  flows,  only  51  grid  points  were 
used  in  the  radial  direction  to  cover  the  radius  of  the  orifice  in  the  first 
problem  and  101  were  used  for  the  second  example. 


6.0  CONCLUDING  REMARKS 

A  random-choice  method  has  been  presented  in  detail  for  obtaining  fairly 
practical  and  efficient  numerical  solutions  for  both  two-dimensional  planar  and 
axisymmetric  steady  supersonic  flows.  Innovative  techniques  were  introduced  for 
solving  the  Riemann  problem  iteratively,  naturally  handling  boundary  conditions 
at  body  and  free-jet  surfaces,  and  computing  only  certain  parts  of  flow  fields 
of  interest.  Many  interesting  and  practical  numerical  solutions  were  also  given 
for  a  variety  of  different  planar-  and  axisymmetric-f low  problems,  and  compared, 
in  most  cases,  with  known  analytical  and  numerical  solutions,  in  order  to  demon¬ 
strate  the  applicability,  capability,  and  limitations  of  this  new  random-choice 
method. 


For  solving  planar-flow  problems  the  random-choice  method  is  limited 
only  in  that  the  flow  must  be  everywhere  supersonic  (see  section  4.1).  Outside 
of  this  single  but  restrictive  limitation  the  method  is  highly  successful  in 
predicting  planar  supersonic  flow  fields  both  efficiently  and  accurately.  For 
the  solution  of  axisymmetric-f low  problems  the  random-choice  method  is  further 
limited  in  that  the  flow  must  be  supersonic  when  solving  the  intermediate  planar 
Riemann  problem  (section  4.2),  before  the  operator-splitting  correction  gives 
the  final  axisymmetric  results.  This  can  be  a  severe  limitation  in  that  flows 
over  most  military  projectiles  cannot  be  solved  by  the  present  random-choice 
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method.  It  is  precisely  these  limitations  that  make  the  random-choice  method 
for  solving  two-dimensional  steady  planar  and  axisymmetr ic  supersonic  flows 
less  extensive  and  therefore  less  useful  in  applications  than  its  sequel  for 
solving  one-dimensional  nonstationary  flows  (for  which  the  latter  has  no  such 
similar  limitations). 

The  problem  of  having  to  resort  to  an  undesirably  fine  grid  mesh  in 
order  to  maintain  numerical  accuracy  when  solving  axisymmetr ic-f low  problems 
that  include  the  axis  of  symmetry  is  the  direct  result  of  inaccuracies  stemming 
from  the  application  of  the  operator-splitting  technique,  as  mentioned  earlier 
in  section  4.2,  and  which  also  became  obvious  from  numerical  results  presented 
in  section  5.2. 

However,  the  time  to  solve  planar  and  axisymmetric  flow  problems  with 
the  random-choice  method  is  not  unreasonable.  Typical  times  to  solve  a  single 
Riemann  problem  are  13  ms  for  the  planar  case  and  15  ms  for  the  axisymmetric 
case.  The  additional  15%  is  due  to  additional  calculations  to  include  the 
operator-splitting  correction.  These  times  are  calculated  by  taking  the  total 
time  to  execute  a  job  divided  by  the  total  number  of  Riemann  problems.  Further, 
these  times  are  given  for  the  UTIAS  computer,  which  is  Perkin-Elmer  32-bit  mini¬ 
computer  with  a  3250  central  processor.  Such  computation  times  would  be  lower 
by  a  factor  of  about  four  or  five  for  a  IBM-3033  or  CDC-6600  machine. 

Future  work  might  involve  the  development  of  a  second-order  random- 
choice  method  for  solving  axisymmetr ic-f low  problems,  one  that  does  not  require 
the  operator-splitting  technique.  In  other  words,  the  Riemann  problem  for  an 
axisymmetric  flow  would  be  solved  directly,  instead  of  the  present  first-order 
indirect  method  of  first  solving  the  planar  Riemann  problem  and  then  utilizing 
the  operator-splitting  technique  to  correct  this  intermediate  result  to  get  the 
ax isymmetric-f low  solution.  The  advantages  of  this  second-order  method  would  be 
two-fold.  Firstly,  the  operator-splitting  technique  would  be  eliminated  along 
with  its  inaccurate  correction  if  the  radius  is  small  or  near  zero.  This  means 
that  a  coarser  mesh  could  be  used  in  the  computations,  although  each  cell  would 
require  more  computation  time.  Secondly,  and  more  importantly,  axisymmetric- 
flow  problems  that  are  unsolvable  by  the  present  first-order  method,  because  a 
solution  to  the  planar  Riemann  problem  did  not  exist,  will  be  able  to  be  solved. 
The  only  limitation  would  be  that  the  flow  must  be  everywhere  supersonic.  Then 
the  second-order  method  would  be  able  to  provide  axisymmetric  supersonic  flow 
solutions  for  most  military  projectiles,  for  example,  but  still  not  over  blunt 
bod  ies  . 
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Fig.  2.  Four  different  wave  patterns  for  the  Riemann  problem, 
composed  of  elemental  oblique  shock  waves  (I?)  and 
Prandtl-Meyer  rarefaction  waves  (^) ,  and  separated  by 
a  slip  stream. 


b) 


Fig.  3.  Diagram  of  leftward  (a)  and  rightward  (b)  rarefaction 

waves  and  the  definition  of  various  flow  and  wave  angles. 


Fig.  5.  Riemann  problem  showing  the  left  state  S-^, 
right  state  Sr,  and  the  continuity  of  pres¬ 
sure  p^  and  flow  angle  across  the  slip 
stream. 
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Fig.  7 


Illustration  of  the  planar  or  axisymmetric  flow  field 
from  a  body  in  a  supersonic  flow,  showing  regions  of 
free-stream  conditions  (A),  computed  flow  conditions 
(B),  and  uncomputed  flow  conditions  (C) . 
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Fig.  8.  Maximum  deflection  angle  of  the  flow  through  an 
oblique  shock  wave  that  produces  a  sonic  flow 
behind  it  (y  =  7/5) . 
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Fig.  9.  Planar  flow  over  a  compressive  corner. 
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Fig.  10a.  Flow  over  an  expansive  corner. 
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Comparison  of  RCM  and  exact  solutions  for  the 
flow  over  an  expansive  corner  (x  =  4.0  units). 
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Fig.  11a.  Flow  over  a  contoured  compressive 
corner  (Mi  =  1.30). 


Fig.  lib.  Flow  over  a  contoured  compressive 
comer  (Mi  =  1.64). 
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Fig.  11c 


Flow  over  a  contoured  compressive 
corner  (Mi  =  2.00). 
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Fig.  lid.  Flow  over  a  contoured  compressive 
corner  (Mj  =  4.00). 
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Flow  in  a  duct  inlet  made  of  a  flat  plate 
and  compressive  wedge  (6  degrees). 
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Flow  in  a  duct  inlet  made  of  two  compressive 
wedges  (2  and  6  degrees). 
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Flow  in  a  duct  inlet  made  of  a  flat  plate 
and  sharp  expansive  corner  (2  degrees). 
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Fig.  16.  Flow  in  a  duct  inlet  made  of  two  compressive 
wedges  (6  and  3  degrees) . 
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Distance  x/L  (L  is  arbitrary) 
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Fig.  20.  Overpressure  signatures  on  the  surface  of,  and  at  1, 
10  and  100  body  thicknesses  from,  a  parabolic  shaped 
airfoil  (t/L  =  0.10)  in  a  supersonic  flow  (Mj  =  2.0). 
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Fig.  21.  Numerical  computations  for  a  free  jet  (Mj  =  1.40). 


RCM  (5) 


T 


Oblique  shock 


Fig.  23.  Illustration  of  supersonic  flow  over 
a  cone  (aligned  with  the  flow) . 
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Fig.  24a.  Comparison  of  RCM  and  exact  results  for  a  supersonic 
flow  over  a  15-degree  cone  (Mj^  =  2.5787,  80  grid 
nodes  downstream  of  the  cone  apex) . 
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Fig.  24b.  Comparison  of  RCM  and  exact  results  for  a  supersonic 
flow  over  a  15-degree  cone  (Mj  =  2.5787,  160  grid 
nodes  downstream  of  the  cone  apex) . 


Fig.  24c.  Comparison  of  RCM  and  exact  results  for  a  supersonic 
flow  over  a  15-degree  cone  (Mj  =  2.5787,  320  grid 
nodes  downstream  of  the  cone  apex) . 
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Fig.  25.  Comparison  of  RCM  and  exact  results  for  the 
pressure  coefficient  for  the  supersonic  flow 
(Mj  =  2.0)  over  a  cone-cylinder. 
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Fig.  26a.  Comparison  of  RCM  and  exact  solutions  for  the 
pressure  coefficient  for  the  supersonic  flow 
(Mj  =  2.0)  over  a  cylindrical  body  with  a 
parabolic  nose. 
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Fig.  26b.  Comparison  of  RCM  and  exact  solutions  for  the 
pressure  coefficient  for  the  supersonic  flow 
(Mj  =  3.0)  over  a  cylindrical  body  with  a 
parabolic  nose. 
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Fig.  26c.  Comparison  of  RCM  and  exact  solutions  for  the 
pressure  coefficient  for  the  supersonic  flow 
(Mi  =  3.92)  over  a  cylindrical  body  with  a 
parabolic  nose. 


ig.  28.  Pressure  profiles  for  a  supersonic  flow  (Mi  =  2.075)  over  a  symmetric 
double-cone  (semivertex  angle  of  10  degrees).  The  RCM  computations 
do  not  include  the  upper  left  and  lower  right  regions. 
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Fig.  30.  Overpressure  signatures  at  radii  of  1,  10,  and 
100  diameters  from  a  parabolic-spindle  shaped 
body  (D/L  =  0.10)  in  a  supersonic  flow  (Mj  = 2) . 
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APPENDIX  A 


COMPUTER-PROGRAM  LISTING 

FOR  SOLVING  TWO-DIMENSIONAL  STEADY  SUPERSONIC  FLOW  PROBLEMS 


In  the  present  numerical  investigation,  many  slightly  different  computer 
programs  were  used  to  obtain  numerical  solutions  for  two-dimensional  planar  and 
axisymmetric  flow  fields  of  an  inviscid,  compressible  and  polytropic  gas  over 
and  around  different  shaped  bodies,  inside  duct  inlets,  and  in  free  jets  from 
slots  and  orifices.  Only  one  typical  version  of  these  computer  programs,  which 
is  based  on  the  random-choice  method  is  listed  here,  that  which  was  used  for 
computations  of  the  flow  over  an  axisymmetric  body. 

The  algorithm  for  the  random-choice  method  solves  a  basic  Riemann  prob¬ 
lem  for  two-dimensional  steady  supersonic  flow  between  adjacent  grid  points  in 
the  radial  direction  at  increments  of  axial  distance  that  are  sufficiently  small 
to  satisfy  the  well-known  Courant-Friedr ichs-Lewy  stability  condition.  It  also 
includes  an  operator-splitting  technique  whi-'h  can  modify  the  initial  Riemann 
solution  for  two-dimensional  planar  flows  so  that  the  desired  solution  for  two- 
dimensional  axisymmetric  flows  can  be  obtained. 

In  the  computer  program,  all  dependent  thermodynamic  and  dynamic  values 
are  nond imensional .  They  are  designated  here  by  using  an  overhead  bar  and  given 
as  follows: 

P  =  P/P0.  P  =  p/p0*  «  =  «/ (p0/Po)1/2'  7  =  v/(p0/p0)1/2, 

where  the  dependent  variables  p,  p,  u,  and  v  denote  the  static  pressure  and 
density  and  the  radial  and  axial  components  of  gas  velocity,  respectively.  The 
variables  with  tho  subscript  o  are  free-stream  or  reference  flow  conditions. 

The  independent  distance  variables  of  x  and  r  are  also  given  in  nondimensional 
form  as  follows: 


x  =  x/L, 


r  =  r/L, 


where  L  is  a  certain  characteristic  reference  length,  which  can  be  a  body  length 
or  diameter. 

The  computer  code  consists  of  one  main  program  followed  by  five  subrou¬ 
tines  and  seven  functions.  The  main  purposes  of  the  subroutines  and  functions 
are  not  described  here  because  they  are  already  outlined  as  commentary  at  the 
beginning  of  each  listing. 
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